lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into compute-mliap

This commit is contained in:
athomps 2020-07-03 14:26:38 -06:00 committed by GitHub
parent 7fe2df423f 558d2eb84f
commit c3f8644613
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
210 changed files with 5111 additions and 3161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.github/CODEOWNERS vendored
View File

@ -119,6 +119,8 @@ src/utils.* @akohlmey @rbberger
tools/msi2lmp/* @akohlmey
tools/emacs/* @HaoZeke
tools/singularity/* @akohlmey @rbberger
tools/code_standard/* @rbberger
tools/valgrind/* @akohlmey
# tests
unittest/* @akohlmey @rbberger

15
cmake/Modules/Packages/GPU.cmake
View File

@ -268,10 +268,17 @@ elseif(GPU_API STREQUAL "HIP")
if(HIP_PLATFORM STREQUAL "hcc")
configure_file(${CU_FILE} ${CU_CPP_FILE} COPYONLY)
add_custom_command(OUTPUT ${CUBIN_FILE}
VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco -t="${HIP_ARCH}" -f=\"-O3 -ffast-math -DUSE_HIP -D_${GPU_PREC_SETTING} -I${LAMMPS_LIB_SOURCE_DIR}/gpu\" -o ${CUBIN_FILE} ${CU_CPP_FILE}
DEPENDS ${CU_CPP_FILE}
COMMENT "Generating ${CU_NAME}.cubin")
if(HIP_COMPILER STREQUAL "clang")
add_custom_command(OUTPUT ${CUBIN_FILE}
VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco --offload-arch=${HIP_ARCH} -O3 -ffast-math -DUSE_HIP -D_${GPU_PREC_SETTING} -I${LAMMPS_LIB_SOURCE_DIR}/gpu -o ${CUBIN_FILE} ${CU_CPP_FILE}
DEPENDS ${CU_CPP_FILE}
COMMENT "Generating ${CU_NAME}.cubin")
else()
add_custom_command(OUTPUT ${CUBIN_FILE}
VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --genco -t="${HIP_ARCH}" -f=\"-O3 -ffast-math -DUSE_HIP -D_${GPU_PREC_SETTING} -I${LAMMPS_LIB_SOURCE_DIR}/gpu\" -o ${CUBIN_FILE} ${CU_CPP_FILE}
DEPENDS ${CU_CPP_FILE}
COMMENT "Generating ${CU_NAME}.cubin")
endif()
elseif(HIP_PLATFORM STREQUAL "nvcc")
add_custom_command(OUTPUT ${CUBIN_FILE}
VERBATIM COMMAND ${HIP_HIPCC_EXECUTABLE} --fatbin --use_fast_math -DUSE_HIP -D_${GPU_PREC_SETTING} ${HIP_CUDA_GENCODE} -I${LAMMPS_LIB_SOURCE_DIR}/gpu -o ${CUBIN_FILE} ${CU_FILE}

2
cmake/Modules/Packages/USER-INTEL.cmake
View File

@ -77,7 +77,7 @@ else()
foreach(_FLAG -O2 -fp-model fast=2 -no-prec-div -qoverride-limits -qopt-zmm-usage=high -qno-offload -fno-alias -ansi-alias -restrict)
check_cxx_compiler_flag("${_FLAG}" COMPILER_SUPPORTS${_FLAG})
if(COMPILER_SUPPORTS${_FLAG})
target_compile_options(lammps PRIVATE ${_FLAG})
target_compile_options(lammps PRIVATE ${_FLAG})
endif()
endforeach()
endif()

2
doc/lammps.1
View File

@ -1,4 +1,4 @@
.TH LAMMPS "15 June 2020" "2020-06-15"
.TH LAMMPS "30 June 2020" "2020-06-30"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator.

2
doc/src/Speed_gpu.rst
View File

@ -88,7 +88,7 @@ GPUs/node to 1.
Using the "-pk" switch explicitly allows for setting of the number of
GPUs/node to use and additional options. Its syntax is the same as
same as the "package gpu" command. See the :doc:`package <package>`
the "package gpu" command. See the :doc:`package <package>`
command doc page for details, including the default values used for
all its options if it is not specified.

10
doc/src/compute_orientorder_atom.rst
View File

@ -61,8 +61,8 @@ The summation is over the *nnn* nearest
neighbors of the central atom.
The angles :math:`theta` and :math:`phi` are the standard spherical polar angles
defining the direction of the bond vector :math:`r_{ij}`.
The phase and sign of :math:`Y_{lm}` follow the standard conventions,
so that :math:`{\rm sign}(Y_{ll}(0,0)) = (-1)^l`.
The phase and sign of :math:`Y_{lm}` follow the standard conventions,
so that :math:`{\rm sign}(Y_{ll}(0,0)) = (-1)^l`.
The second equation defines :math:`Q_l`, which is a
rotationally invariant non-negative amplitude obtained by summing
over all the components of degree *l*\ .
@ -181,13 +181,13 @@ values for each atom will be added to the output array, which are real numbers.
If the keyword *components* is set, then the real and imaginary parts
of each component of *normalized* :math:`\hat{Y}_{lm}` will be added to the
output array in the following order: :math:`{\rm Re}(\hat{Y}_{-m}), {\rm Im}(\hat{Y}_{-m}),
{\rm Re}(\hat{Y}_{-m+1}), {\rm Im}(\hat{Y}_{-m+1}), \dots , {\rm Re}(\hat{Y}_m), {\rm Im}(\hat{Y}_m)`.
output array in the following order: :math:`{\rm Re}(\hat{Y}_{-m}), {\rm Im}(\hat{Y}_{-m}),
{\rm Re}(\hat{Y}_{-m+1}), {\rm Im}(\hat{Y}_{-m+1}), \dots , {\rm Re}(\hat{Y}_m), {\rm Im}(\hat{Y}_m)`.
In summary, the per-atom array will contain *nlvalues* columns, followed by
an additional *nlvalues* columns if *wl* is set to yes, followed by
an additional *nlvalues* columns if *wl/hat* is set to yes, followed
by an additional 2\*(2\* *ldegree*\ +1) columns if the *components*
by an additional 2\*(2\* *ldegree*\ +1) columns if the *components*
keyword is set.
These values can be accessed by any command that uses per-atom values

14
doc/src/pair_coeff.rst
View File

@ -89,14 +89,20 @@ LAMMPS so long as they are in the format LAMMPS expects, as discussed
on the individual doc pages. The first line of potential files may
contain metadata with upper case tags followed their value. These may
be parsed and used by LAMMPS. Currently supported are the "DATE:"
tag and the "UNITS:" tag. For pair styles that have been programmed
tag and the ``UNITS:`` tag. For pair styles that have been programmed
to support the metadata, the value of the "DATE:" tag is printed to
the screen and logfile so that the version of a potential file can be
later identified. The "UNITS:" tag indicates the :doc:`units <units>`
later identified. The ``UNITS:`` tag indicates the :doc:`units <units>`
setting required for this particular potential file. If the potential
file ware created for a different sets of units, LAMMPS will terminate
file was created for a different sets of units, LAMMPS will terminate
with an error. If the potential file does not contain the tag, no
check will be made.
check will be made and it is the responsibility of the user to determine
that the unit style is correct.
In some select cases and for specific combinations of unit styles,
LAMMPS is capable of automatically converting potential parameters
from a file. In those cases, a warning message signaling that an
automatic conversion has happened is printed to the screen.
When a pair_coeff command using a potential file is specified, LAMMPS
looks for the potential file in 2 places. First it looks in the

2
doc/src/pair_eam.rst
View File

@ -76,7 +76,7 @@ Examples
pair_style eam
pair_coeff * * cuu3
pair_coeff 1*3 1\*3 niu3.eam
pair_coeff 1*3 1*3 niu3.eam
pair_style eam/alloy
pair_coeff * * ../potentials/NiAlH_jea.eam.alloy Ni Al Ni Ni

4
doc/src/pair_eim.rst
View File

@ -48,7 +48,7 @@ and :math:`sigma_i` are calculated as
\sigma_i = & \sum_{j=i_1}^{i_N} q_j \cdot \psi_{ij} \left(r_{ij}\right) \\
E_i\left(q_i,\sigma_i\right) = & \frac{1}{2} \cdot q_i \cdot \sigma_i
where :math:`\eta_{ji} is a pairwise function describing electron flow from atom
where :math:`\eta_{ji}` is a pairwise function describing electron flow from atom
I to atom J, and :math:`\psi_{ij}` is another pairwise function. The multi-body
nature of the EIM potential is a result of the embedding energy term.
A complete list of all the pair functions used in EIM is summarized
@ -63,7 +63,7 @@ below
\right.\\
\eta_{ji} = & A_{\eta,ij}\left(\chi_j-\chi_i\right)f_c\left(r,r_{s,\eta,ij},r_{c,\eta,ij}\right) \\
\psi_{ij}\left(r\right) = & A_{\psi,ij}\exp\left(-\zeta_{ij}r\right)f_c\left(r,r_{s,\psi,ij},r_{c,\psi,ij}\right) \\
f_{c}\left(r,r_p,r_c\right) = & 0.510204 \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204
f_{c}\left(r,r_p,r_c\right) = & 0.510204 \cdot \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204
Here :math:`E_b, r_e, r_(c,\phi), \alpha, \beta, A_(\psi), \zeta, r_(s,\psi),
r_(c,\psi), A_(\eta), r_(s,\eta), r_(c,\eta), \chi,` and pair function type

30
doc/src/pair_mliap.rst
View File

@ -28,7 +28,7 @@ Examples
.. code-block:: LAMMPS
pair_style mliap model linear InP.mliap.model descriptor sna InP.mliap.descriptor
pair_style mliap model quadratic W.mliap.model descriptor sna W.mliap.descriptor
pair_style mliap model quadratic W.mliap.model descriptor sna W.mliap.descriptor
pair_coeff * * In P
Description
@ -38,9 +38,9 @@ Pair style *mliap* provides a general interface to families of
machine-learning interatomic potentials. It allows separate
definitions of the interatomic potential functional form (*model*)
and the geometric quantities that characterize the atomic positions
(*descriptor*). By defining *model* and *descriptor* separately,
(*descriptor*). By defining *model* and *descriptor* separately,
it is possible to use many different models with a given descriptor,
or many different descriptors with a given model. Currently, the pair_style
or many different descriptors with a given model. Currently, the pair_style
supports just two models, *linear* and *quadratic*,
and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap <pair_snap>`, including the linear, quadratic,
and chem variants. Work is currently underway to extend
@ -52,11 +52,11 @@ can be accessed using the related :doc:`compute mliap <compute_mliap>` command.
The pair_style *mliap* command must be followed by two keywords
*model* and *descriptor* in either order. A single
*pair_coeff* command is also required. The first 2 arguments
*pair_coeff* command is also required. The first 2 arguments
must be \* \* so as to span all LAMMPS atom types.
This is followed by a list of N arguments
that specify the mapping of MLIAP
element names to LAMMPS atom types,
element names to LAMMPS atom types,
where N is the number of LAMMPS atom types.
The *model* keyword is followed by a model style, currently limited to
@ -68,7 +68,7 @@ It may contain parameters for many elements. The only requirement is that it
contain at least those element names appearing in the
*pair_coeff* command.
The top of the model file can contain any number of blank and comment lines (start with #),
The top of the model file can contain any number of blank and comment lines (start with #),
but follows a strict format after that. The first non-blank non-comment
line must contain two integers:
@ -79,20 +79,20 @@ This is followed by one block for each of the *nelem* elements.
Each block consists of *nparams* parameters, one per line.
Note that this format is similar, but not identical to that used
for the :doc:`pair_style snap <pair_snap>` coefficient file.
Specifically, the line containing the element weight and radius is omitted,
Specifically, the line containing the element weight and radius is omitted,
since these are handled by the *descriptor*.
The *descriptor* keyword is followed by a descriptor style, and additional arguments.
Currently the only descriptor style is *sna*, indicating the bispectrum component
descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of
Currently the only descriptor style is *sna*, indicating the bispectrum component
descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of
:doc:`pair_style snap <pair_snap>`.
The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped
out by the LAMMPS command parser. A single additional argument specifies the descriptor filename
containing the parameters and setting used by the SNAP descriptor.
The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped
out by the LAMMPS command parser. A single additional argument specifies the descriptor filename
containing the parameters and setting used by the SNAP descriptor.
The descriptor filename usually ends in the *.mliap.descriptor* extension.
The SNAP descriptor file closely follows the format of the
:doc:`pair_style snap <pair_snap>` parameter file.
The SNAP descriptor file closely follows the format of the
:doc:`pair_style snap <pair_snap>` parameter file.
The file can contain blank and comment lines (start
with #) anywhere. Each non-blank non-comment line must contain one
keyword/value pair. The required keywords are *rcutfac* and
@ -102,7 +102,7 @@ In addition, the SNAP descriptor file must contain
the *nelems*, *elems*, *radelems*, and *welems* keywords.
The *nelems* keyword specifies the number of elements
provided in the other three keywords.
The *elems* keyword is followed by a list of *nelems*
The *elems* keyword is followed by a list of *nelems*
element names that must include the element
names appearing in the *pair_coeff* command,
but can contain other names too.

28
doc/src/pair_snap.rst
View File

@ -25,9 +25,9 @@ Description
"""""""""""
Pair style *snap* defines the spectral
neighbor analysis potential (SNAP), a machine-learning
neighbor analysis potential (SNAP), a machine-learning
interatomic potential :ref:`(Thompson) <Thompson20142>`.
Like the GAP framework of Bartok et al. :ref:`(Bartok2010) <Bartok20102>`,
Like the GAP framework of Bartok et al. :ref:`(Bartok2010) <Bartok20102>`,
SNAP uses bispectrum components
to characterize the local neighborhood of each atom
in a very general way. The mathematical definition of the
@ -139,7 +139,7 @@ The SNAP parameter file can contain blank and comment lines (start
with #) anywhere. Each non-blank non-comment line must contain one
keyword/value pair. The required keywords are *rcutfac* and
*twojmax*\ . Optional keywords are *rfac0*\ , *rmin0*\ ,
*switchflag*\ , *bzeroflag*\ , *quadraticflag*\ , *chemflag*\ ,
*switchflag*\ , *bzeroflag*\ , *quadraticflag*\ , *chemflag*\ ,
*bnormflag*\ , *wselfallflag*\ , and *chunksize*\ .
The default values for these keywords are
@ -154,34 +154,34 @@ The default values for these keywords are
* *wselfallflag* = 0
* *chunksize* = 2000
If *quadraticflag* is set to 1, then the SNAP energy expression includes additional quadratic terms
If *quadraticflag* is set to 1, then the SNAP energy expression includes additional quadratic terms
that have been shown to increase the overall accuracy of the potential without much increase
in computational cost :ref:`(Wood) <Wood20182>`.
in computational cost :ref:`(Wood) <Wood20182>`.
.. math::
E^i_{SNAP}(\mathbf{B}^i) = \beta^{\mu_i}_0 + \boldsymbol{\beta}^{\mu_i} \cdot \mathbf{B}_i + \frac{1}{2}\mathbf{B}^t_i \cdot \boldsymbol{\alpha}^{\mu_i} \cdot \mathbf{B}_i
where :math:`\mathbf{B}_i` is the *K*-vector of bispectrum components,
:math:`\boldsymbol{\beta}^{\mu_i}` is the *K*-vector of linear coefficients
for element :math:`\mu_i`, and :math:`\boldsymbol{\alpha}^{\mu_i}`
where :math:`\mathbf{B}_i` is the *K*-vector of bispectrum components,
:math:`\boldsymbol{\beta}^{\mu_i}` is the *K*-vector of linear coefficients
for element :math:`\mu_i`, and :math:`\boldsymbol{\alpha}^{\mu_i}`
is the symmetric *K* by *K* matrix of quadratic coefficients.
The SNAP element file should contain *K*\ (\ *K*\ +1)/2 additional coefficients
for each element, the upper-triangular elements of :math:`\boldsymbol{\alpha}^{\mu_i}`.
If *chemflag* is set to 1, then the energy expression is written in terms of explicit multi-element bispectrum
components indexed on ordered triplets of elements, which has been shown to increase the ability of the SNAP
potential to capture energy differences in chemically complex systems,
potential to capture energy differences in chemically complex systems,
at the expense of a significant increase in computational cost :ref:`(Cusentino) <Cusentino20202>`.
.. math::
E^i_{SNAP}(\mathbf{B}^i) = \beta^{\mu_i}_0 + \sum_{\kappa,\lambda,\mu} \boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i} \cdot \mathbf{B}^{\kappa\lambda\mu}_i
E^i_{SNAP}(\mathbf{B}^i) = \beta^{\mu_i}_0 + \sum_{\kappa,\lambda,\mu} \boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i} \cdot \mathbf{B}^{\kappa\lambda\mu}_i
where :math:`\mathbf{B}^{\kappa\lambda\mu}_i` is the *K*-vector of bispectrum components
for neighbors of elements :math:`\kappa`, :math:`\lambda`, and :math:`\mu` and
:math:`\boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i}` is the corresponding *K*-vector
of linear coefficients for element :math:`\mu_i`. The SNAP element file should contain
where :math:`\mathbf{B}^{\kappa\lambda\mu}_i` is the *K*-vector of bispectrum components
for neighbors of elements :math:`\kappa`, :math:`\lambda`, and :math:`\mu` and
:math:`\boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i}` is the corresponding *K*-vector
of linear coefficients for element :math:`\mu_i`. The SNAP element file should contain
a total of :math:`K N_{elem}^3` coefficients for each of the :math:`N_{elem}` elements.
The keyword *chunksize* is only applicable when using the

2
doc/src/pair_tersoff_mod.rst
View File

@ -74,7 +74,7 @@ formulation of the V_ij term, where it contains an additional c0 term.
.. math::
V_{ij} & = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right]
V_{ij} = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right]
The modified cutoff function :math:`f_C` proposed by :ref:`(Murty) <Murty>` and
having a continuous second-order differential is employed. The

9
doc/src/units.rst
View File

@ -49,6 +49,15 @@ new units. And you must correctly convert all output from the new
units to the old units when comparing to the original results. That
is often not simple to do.
Potential or table files may have a ``UNITS:`` tag included in the
first line indicating the unit style those files were created for.
If the tag exists, its value will be compared to the chosen unit style
and LAMMPS will stop with an error message if there is a mismatch.
In some select cases and for specific combinations of unit styles,
LAMMPS is capable of automatically converting potential parameters
from a file. In those cases, a warning message signaling that an
automatic conversion has happened is printed to the screen.
----------
For style *lj*\ , all quantities are unitless. Without loss of

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -1760,6 +1760,7 @@ mem
memalign
MEMALIGN
membered
memcheck
Mendelev
mer
Meremianin

1
examples/README
View File

@ -88,6 +88,7 @@ melt: rapid melt of 3d LJ system
message: client/server coupling of 2 codes
micelle: self-assembly of small lipid-like molecules into 2d bilayers
min: energy minimization of 2d LJ melt
mliap: examples for using several bundled MLIAP potentials
mscg: parameterize a multi-scale coarse-graining (MSCG) model
msst: MSST shock dynamics
nb3b: use of nonbonded 3-body harmonic pair style

20
examples/USER/qtb/alpha_quartz_qbmsst/alpha_quartz_potential.mod
View File

@ -1,15 +1,15 @@
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
kspace_style pppm 1.0e-4
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000

55
examples/USER/qtb/alpha_quartz_qbmsst/alpha_quartz_qbmsst.in
View File

@ -1,55 +0,0 @@
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 100
timestep ${delta_t}
run 1000
write_restart restart.1000
#Compression restart 2
#Read restart file and load potential again
clear
read_restart restart.1000
include alpha_quartz_potential.mod
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 100
timestep ${delta_t}
restart 1000 restart
run 10000 #10 ps

90
examples/USER/qtb/alpha_quartz_qbmsst/alpha_quartz_qtb.mod
View File

@ -3,60 +3,60 @@
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
include alpha_quartz_potential.mod
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 100
run 2000 # 2 ps
unfix quartz_qtb
unfix scapegoat_qtb
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
unfix quartz_qtb
unfix scapegoat_qtb

55
examples/USER/qtb/alpha_quartz_qbmsst/in.alpha_quartz_qbmsst Normal file
View File

@ -0,0 +1,55 @@
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
run 1000
write_restart restart.1000
#Compression restart 2
#Read restart file and load potential again
clear
read_restart restart.1000
include alpha_quartz_potential.mod
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
restart 1000 restart
run 10000 #10 ps

430
examples/USER/qtb/alpha_quartz_qbmsst/log.15Jun20.alpha_quartz_qbmsst.g++.1 Normal file
View File

@ -0,0 +1,430 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 1 by 1 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 1 ${z_rep}
replicate 2 1 4
orthogonal box = (0.0 0.0 0.0) to (9.832 8.514762 21.6216)
1 by 1 by 1 MPI processor grid
144 atoms
replicate CPU = 0.000271082 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:471)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.301598
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.00117056
estimated relative force accuracy = 8.12908e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 80.09 | 80.09 | 80.09 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34026.791 -2793.6042 1810.0985 9.832 8.514762 21.6216 -37470.578 -37470.432 -27139.363 -6.4345368e-12 0.94245783 4.2212262e-10
200 170.7381 43248.332 -2790.8398 1879.164 9.9554912 8.6217086 21.89317 39337.624 42979.126 47428.246 324.91326 454.85872 -2034.6053
400 258.09921 -28257.8 -2788.3487 1856.1432 9.9146707 8.5863569 21.803402 -19478.873 -29571.375 -35723.152 4633.9026 8487.8103 -626.12005
600 277.77032 -22751.351 -2786.2715 1866.9783 9.9339253 8.6030319 21.845744 -21727.335 -29200.027 -17326.692 -4327.8571 -8218.4965 252.30681
800 349.8665 30508.003 -2784.2204 1873.4953 9.9454706 8.6130304 21.871134 29929.055 33562.672 28032.281 -3188.5605 12329.482 7558.5678
1000 373.67651 -18839.569 -2783.2178 1855.5937 9.9136922 8.5855095 21.80125 -18063.486 -22288.321 -16166.902 -416.09547 -10368.975 9030.4208
1200 423.3474 6846.9905 -2781.9271 1896.2131 9.9855083 8.6477041 21.959181 2147.3938 11765.857 6627.7202 -7627.6782 -1297.6517 -4758.4746
1400 418.54527 -6416.7506 -2781.4358 1834.2719 9.8755745 8.5524986 21.717425 5693.0543 -19487.901 -5455.405 827.66513 -523.1508 -3890.9919
1600 429.42796 3939.8836 -2780.5861 1895.8859 9.984934 8.6472068 21.957918 3755.6959 -1326.4343 9390.3893 1948.1153 4489.8629 1466.0914
1800 447.7623 -8344.6306 -2780.1071 1858.4925 9.9188518 8.5899779 21.812596 -17549.498 3336.8135 -10821.208 1643.4226 -644.56065 -8935.9666
2000 438.1306 -6691.4691 -2780.7407 1871.3547 9.9416812 8.6097487 21.862801 -6959.2196 -8486.8466 -4628.341 -1019.9006 443.03694 -2751.917
Loop time of 2.46815 on 1 procs for 2000 steps with 144 atoms
Performance: 70.012 ns/day, 0.343 hours/ns, 810.323 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.0003 | 2.0003 | 2.0003 | 0.0 | 81.04
Kspace | 0.20006 | 0.20006 | 0.20006 | 0.0 | 8.11
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.023753 | 0.023753 | 0.023753 | 0.0 | 0.96
Output | 0.0001986 | 0.0001986 | 0.0001986 | 0.0 | 0.01
Modify | 0.23896 | 0.23896 | 0.23896 | 0.0 | 9.68
Other | | 0.004907 | | | 0.20
Nlocal: 144 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3943 ave 3943 max 3943 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 41952 ave 41952 max 41952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
timestep 0.001
run 1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.303132
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.00104699
estimated relative force accuracy = 7.27093e-05
using double precision FFTW3
3d grid and FFT values/proc = 5520 1152
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Fix QBMSST v0 = 1.87135e+03
Fix QBMSST p0 = -4.62948e+03
Fix QBMSST e0 = to be -2.78074e+03
Fix QBMSST initial strain rate of -4.01096e-01 established by reducing temperature by factor of 5.00000e-02
Per MPI rank memory allocation (min/avg/max) = 80.1 | 80.1 | 80.1 Mbytes
Step v_T_qm Press TotEng Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
0 300 -6922.9433 -2780.7394 1871.3547 9.9416812 8.6097487 21.862801 -4819.9907 10.953265 -190.51273 0 0
200 294.95797 54876.416 -2779.2988 1723.7621 9.9416812 8.6097487 20.138495 108897.19 -29.773973 -9271.7281 6.1518102 -15.057867
400 288.3711 139521.03 -2778.7321 1628.5574 9.9416812 8.6097487 19.026231 222107.71 8.0682073 24727.575 10.120041 -28.714693
600 280.56521 98070.281 -2779.8934 1687.2434 9.9416812 8.6097487 19.711852 164558.51 2.6076928 16005.656 7.6739491 -42.705007
800 274.94701 106060.26 -2779.2916 1651.0723 9.9416812 8.6097487 19.289269 176842.6 -39.645354 -1804.9466 9.1815975 -56.628078
1000 268.47106 189695.34 -2779.4951 1492.6355 9.9416812 8.6097487 17.438272 277351.5 -84.834482 -33116.996 15.785409 -69.870519
Loop time of 2.05219 on 1 procs for 1000 steps with 144 atoms
Performance: 42.101 ns/day, 0.570 hours/ns, 487.284 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.6815 | 1.6815 | 1.6815 | 0.0 | 81.94
Kspace | 0.10373 | 0.10373 | 0.10373 | 0.0 | 5.05
Neigh | 0.0061183 | 0.0061183 | 0.0061183 | 0.0 | 0.30
Comm | 0.012444 | 0.012444 | 0.012444 | 0.0 | 0.61
Output | 0.00014687 | 0.00014687 | 0.00014687 | 0.0 | 0.01
Modify | 0.24529 | 0.24529 | 0.24529 | 0.0 | 11.95
Other | | 0.002948 | | | 0.14
Nlocal: 144 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4243 ave 4243 max 4243 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 48210 ave 48210 max 48210 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 48210
Ave neighs/atom = 334.792
Neighbor list builds = 8
Dangerous builds = 0
write_restart restart.1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.306435
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.000955688
estimated relative force accuracy = 6.63689e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 3
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
#Compression restart 2
#Read restart file and load potential again
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.1000
restoring atom style charge from restart
orthogonal box = (-0.05484062286382799 -0.04749337384227555 2.0916641327653274) to (9.886840622863804 8.562255373842252 19.52993586723476)
1 by 1 by 1 MPI processor grid
restoring pair style hybrid/overlay from restart
144 atoms
read_restart CPU = 0.0002563 secs
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:471)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
Resetting global fix info from restart file:
fix style: qbmsst, fix ID: shock
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
timestep 0.001
restart 1000 restart
run 10000 #10 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.306435
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.000955688
estimated relative force accuracy = 6.63689e-05
using double precision FFTW3
3d grid and FFT values/proc = 5280 1080
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 3
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 80.12 | 80.12 | 80.12 Mbytes
Step v_T_qm Press TotEng Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
1000 268.47106 189686.77 -2781.5194 1492.6355 9.9416812 8.6097487 17.438272 277378.37 -84.692548 -33090.129 15.785409 0
1500 362.13476 692245.96 -2800.9352 1011.2037 9.9416812 8.6097487 11.813766 661095.53 188.71833 -49928.712 35.851981 -24.11484
2000 860.78914 714816.8 -2830.893 997.64749 9.9416812 8.6097487 11.65539 653537.64 852.68158 -68765.537 36.41702 -44.978484
2500 1620.8281 709511.19 -2840.8217 1000.3425 9.9416812 8.6097487 11.686875 660030.01 1184.3105 -60030.892 36.304689 -65.69966
3000 2395.6824 649526.84 -2832.6859 995.56591 9.9416812 8.6097487 11.631071 660984.37 939.07209 -63050.693 36.503782 -86.383242
3500 3034.6774 715794.56 -2822.6098 995.8622 9.9416812 8.6097487 11.634532 712849.74 1055.7295 -10938.816 36.491433 -106.99315
4000 3487.9039 736791.25 -2804.1216 994.13867 9.9416812 8.6097487 11.614397 765817.85 943.15747 40595.305 36.563271 -127.76315
4500 3718.6279 813775.8 -2788.1942 995.82514 9.9416812 8.6097487 11.634099 881961.06 1370.5559 158141.68 36.492977 -148.68649
5000 3691.4947 750146.58 -2770.5541 1018.4785 9.9416812 8.6097487 11.898756 770500.36 196.2793 65528.786 35.548762 -169.8589
5500 3585.8602 831522.51 -2766.0198 1005.6834 9.9416812 8.6097487 11.749273 916093.67 1088.1987 200476.48 36.082073 -190.89436
6000 3431.6405 749891.94 -2771.6404 1011.9077 9.9416812 8.6097487 11.82199 781321.11 268.24344 70882.55 35.82264 -212.20913
6500 3350.2876 666113.16 -2780.4124 1028.8353 9.9416812 8.6097487 12.019753 749294.32 371.38231 52939.676 35.117081 -233.59556
7000 3339.2397 675783.2 -2782.7559 1022.6541 9.9416812 8.6097487 11.947539 690109.39 -26.949124 -11388.054 35.374719 -254.95868
7500 3395.582 726601.74 -2784.7652 1018.1439 9.9416812 8.6097487 11.894846 759167.86 506.5811 53917.852 35.56271 -276.24361
8000 3393.2372 625141.93 -2771.6398 1035.4915 9.9416812 8.6097487 12.097517 598674.46 -895.80046 -92142.112 34.839641 -297.61681
8500 3272.9752 659367.77 -2776.608 1031.8188 9.9416812 8.6097487 12.054609 688358.42 -142.30814 -5513.8593 34.992722 -318.94541
9000 3277.8848 724828.76 -2777.6502 1017.6314 9.9416812 8.6097487 11.888859 724452.11 58.574942 18775.738 35.58407 -340.1718
9500 3273.7854 620652.38 -2780.0794 1023.5922 9.9416812 8.6097487 11.958499 747175.42 317.3826 46458.505 35.335617 -361.41643
10000 3329.1766 668606.38 -2786.3493 1022.9534 9.9416812 8.6097487 11.951035 703351.81 168.14538 2103.38 35.362244 -382.64609
10500 3398.9956 642919.16 -2784.2833 1016.2587 9.9416812 8.6097487 11.872822 661298.16 -230.03577 -45520.34 35.641287 -403.78721
11000 3418.7053 675754.06 -2782.6318 1005.7483 9.9416812 8.6097487 11.75003 689789.84 -136.97148 -25773.422 36.079372 -424.97556
Loop time of 32.4277 on 1 procs for 10000 steps with 144 atoms
Performance: 26.644 ns/day, 0.901 hours/ns, 308.378 timesteps/s
99.5% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 28.397 | 28.397 | 28.397 | 0.0 | 87.57
Kspace | 1.0225 | 1.0225 | 1.0225 | 0.0 | 3.15
Neigh | 0.27594 | 0.27594 | 0.27594 | 0.0 | 0.85
Comm | 0.1797 | 0.1797 | 0.1797 | 0.0 | 0.55
Output | 0.10409 | 0.10409 | 0.10409 | 0.0 | 0.32
Modify | 2.4112 | 2.4112 | 2.4112 | 0.0 | 7.44
Other | | 0.03707 | | | 0.11
Nlocal: 144 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 5541 ave 5541 max 5541 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 74662 ave 74662 max 74662 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 74662
Ave neighs/atom = 518.486
Neighbor list builds = 207
Dangerous builds = 0
Total wall time: 0:00:37

430
examples/USER/qtb/alpha_quartz_qbmsst/log.15Jun20.alpha_quartz_qbmsst.g++.4 Normal file
View File

@ -0,0 +1,430 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first uses fix qtb to equilibrate alpha quartz structure to an initial state with quantum nuclear correction and then simulate shock induced phase transition through the quantum thermal bath multi-scale shock technique
variable x_rep equal 2 #plot is made with x_rep = 8 #x-direction replication number
variable y_rep equal 1 #plot is made with y_rep = 5 #y-direction replication number
variable z_rep equal 4 #plot is made with z_rep = 15 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
variable v_msst equal 78.0 #Shock velocity (Angstrom/ps in metal units)
variable q_msst equal 40.0 #Box mass-like parameter in the MSST (mass^2/length^4, where mass=grams/mole and length=Angstrom in metal units)
variable tscale_msst equal 0.05 #Temperature reduction parameter in the MSST (unitless)
variable eta_qbmsst equal 1.0 #Coupling constant between the shock and the quantum thermal bath (unitless constant)
## This part uses fix qtb to prepare alpha-quartz with quantum nuclear correction of the initial state
include alpha_quartz_qtb.mod
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 2 by 2 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 1 ${z_rep}
replicate 2 1 4
orthogonal box = (0.0 0.0 0.0) to (9.832 8.514762 21.6216)
1 by 1 by 4 MPI processor grid
144 atoms
replicate CPU = 0.000225782 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:471)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 200
run 2000 # 2 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.301598
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.00117056
estimated relative force accuracy = 8.12908e-05
using double precision FFTW3
3d grid and FFT values/proc = 2400 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 79.7 | 79.7 | 79.7 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34026.791 -2793.6042 1810.0985 9.832 8.514762 21.6216 -37470.578 -37470.432 -27139.363 1.0530512e-10 0.94245783 4.0087238e-10
200 153.57631 45538.205 -2790.8177 1873.0866 9.9447472 8.612404 21.869543 41721.016 44095.248 50798.351 -3961.4596 1223.325 2871.656
400 234.74785 -34404.175 -2789.0189 1850.2127 9.9041 8.5772024 21.780156 -28329.333 -39376.313 -35506.88 -1154.5043 -5411.1071 2246.6749
600 265.24833 -20905.145 -2786.2727 1874.9981 9.948129 8.6153326 21.87698 -22753.886 -21091.083 -18870.468 -4645.5548 2968.2945 1415.0311
800 297.79035 32990.58 -2784.8247 1853.6946 9.910309 8.5825796 21.79381 30061.364 35359.18 33551.195 -3092.2971 1525.52 -6461.0249
1000 367.71884 -27539.239 -2783.0102 1864.7161 9.9299114 8.5995557 21.836917 -20273.387 -38720.429 -23623.901 7639.0334 -866.35665 543.52723
1200 399.77109 3807.7814 -2781.511 1893.4978 9.9807399 8.6435745 21.948695 1625.8226 7441.2236 2356.298 -4057.1674 3814.9305 1528.4567
1400 466.57962 -4148.235 -2780.1546 1851.5925 9.9065614 8.5793341 21.785568 -10883.19 1816.768 -3378.2828 896.25296 -7208.541 -42.253127
1600 497.86539 14505.31 -2778.9409 1882.2616 9.9609584 8.6264432 21.905193 8268.1103 20614.738 14633.082 -2690.5669 6807.3187 11995.878
1800 557.31182 -108.04462 -2778.1875 1875.514 9.9490413 8.6161228 21.878986 948.68308 -1929.7575 656.94053 -1628.2172 -6594.5909 -4423.4368
2000 480.39449 -8852.2243 -2778.4963 1862.9552 9.9267847 8.596848 21.830042 -18274.307 3038.8369 -11321.203 -5002.1016 12023.282 6845.2769
Loop time of 1.42181 on 4 procs for 2000 steps with 144 atoms
Performance: 121.535 ns/day, 0.197 hours/ns, 1406.656 timesteps/s
87.5% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.72578 | 0.80093 | 0.87518 | 6.1 | 56.33
Kspace | 0.33737 | 0.41245 | 0.48642 | 8.4 | 29.01
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.066098 | 0.071334 | 0.076039 | 1.6 | 5.02
Output | 0.00021172 | 0.00039291 | 0.00093484 | 0.0 | 0.03
Modify | 0.090105 | 0.1077 | 0.11384 | 3.1 | 7.58
Other | | 0.029 | | | 2.04
Nlocal: 36 ave 36 max 36 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2614 ave 2614 max 2614 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 10488 ave 11326 max 9404 min
Histogram: 1 0 0 0 0 0 2 0 0 1
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
## This part demonstrates how to retart fix qbmsst during any stage of the shock simulation.
## PPPM may break down when compression ratio becomes extremely large. One can always use this restart technique to resume the shock simulation.
#Compression restart 1
reset_timestep 0
#Beta is the number of time steps between each update of the quantum bath temperature. Setting a larger beta can reduce thermal flactuations.
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale ${tscale_msst} damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.05 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 200
timestep ${delta_t}
timestep 0.001
run 1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.30088
grid = 9 8 15
stencil order = 5
estimated absolute RMS force accuracy = 0.00120534
estimated relative force accuracy = 8.37062e-05
using double precision FFTW3
3d grid and FFT values/proc = 2400 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Fix QBMSST v0 = 1.86296e+03
Fix QBMSST p0 = -1.13219e+04
Fix QBMSST e0 = to be -2.77850e+03
Fix QBMSST initial strain rate of -4.21890e-01 established by reducing temperature by factor of 5.00000e-02
Per MPI rank memory allocation (min/avg/max) = 79.7 | 79.7 | 79.7 Mbytes
Step v_T_qm Press TotEng Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
0 300 -9106.318 -2778.4963 1862.9552 9.9267847 8.596848 21.830042 -11562.002 12.009862 -240.0699 0 0
200 296.47213 25984.111 -2777.5178 1770.2164 9.9267847 8.596848 20.743332 64970.204 -25.305765 -1564.7673 3.8828772 -15.16768
400 291.06707 69977.517 -2777.6325 1684.893 9.9267847 8.596848 19.743515 144833.82 -12.184734 6667.384 7.4552796 -29.607028
600 287.21118 39706.699 -2778.0322 1716.9533 9.9267847 8.596848 20.119196 87971.152 -38.593844 -23279.741 6.1129484 -43.751298
800 284.33611 18833.281 -2778.1637 1792.7576 9.9267847 8.596848 21.007468 43725.433 -8.1267799 -3885.5802 2.9391022 -58.454556
1000 281.98328 -6030.6935 -2778.3314 1881.8369 9.9267847 8.596848 22.051297 -14118.602 1.3183874 13055.078 -0.79055793 -73.780965
Loop time of 1.25215 on 4 procs for 1000 steps with 144 atoms
Performance: 69.001 ns/day, 0.348 hours/ns, 798.628 timesteps/s
90.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.67979 | 0.73665 | 0.8091 | 5.4 | 58.83
Kspace | 0.18687 | 0.25893 | 0.31544 | 9.1 | 20.68
Neigh | 0.0011306 | 0.0012404 | 0.0013735 | 0.3 | 0.10
Comm | 0.040339 | 0.041345 | 0.042296 | 0.4 | 3.30
Output | 0.00020051 | 0.00035506 | 0.00081801 | 0.0 | 0.03
Modify | 0.19595 | 0.2007 | 0.20253 | 0.6 | 16.03
Other | | 0.01292 | | | 1.03
Nlocal: 36 ave 38 max 34 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 2527.75 ave 2547 max 2518 min
Histogram: 2 0 0 1 0 0 0 0 0 1
Neighs: 10194.8 ave 11177 max 9437 min
Histogram: 2 0 0 0 0 0 1 0 0 1
Total # of neighbors = 40779
Ave neighs/atom = 283.188
Neighbor list builds = 6
Dangerous builds = 0
write_restart restart.1000
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.302953
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.00105569
estimated relative force accuracy = 7.33134e-05
using double precision FFTW3
3d grid and FFT values/proc = 2640 288
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
#Compression restart 2
#Read restart file and load potential again
clear
using 1 OpenMP thread(s) per MPI task
read_restart restart.1000
restoring atom style charge from restart
orthogonal box = (-0.04739235907204603 -0.041042988010289584 -0.21484841641189512) to (9.879392359072014 8.555804988010294 21.83644841641206)
1 by 1 by 4 MPI processor grid
restoring pair style hybrid/overlay from restart
144 atoms
read_restart CPU = 0.000472307 secs
include alpha_quartz_potential.mod
#This script implements the BKS pair potential for various silicon dioxide compounds. Inner part is fixed with a harmonic potential. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style hybrid/overlay buck/coul/long ${cut_off} table linear 39901
pair_style hybrid/overlay buck/coul/long 10 table linear 39901
pair_coeff 1 1 buck/coul/long 0.0 1.0 0.0 #No interactions between Si atoms
pair_coeff 1 2 buck/coul/long 18003.757200 0.205205 133.538100
pair_coeff 2 2 buck/coul/long 1388.773000 0.362319 175.000000 #BKS interaction in PRL 64 1955 (1990)
pair_modify shift yes
pair_coeff 1 2 table potential_SiO2.TPF Si-O ${cut_off}
pair_coeff 1 2 table potential_SiO2.TPF Si-O 10
pair_coeff 2 2 table potential_SiO2.TPF O-O ${cut_off} #See the potential file for more information
pair_coeff 2 2 table potential_SiO2.TPF O-O 10
WARNING: 1 of 39901 force values in table are inconsistent with -dE/dr.
Should only be flagged at inflection points (src/pair_table.cpp:471)
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
#Use the same fix id and add no tscale if the system is already compressed
fix shock all qbmsst z ${v_msst} q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q ${q_msst} tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp ${damp_qtb} f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta ${eta_qbmsst} beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init ${temperature}
fix shock all qbmsst z 78 q 40 tscale 0.0 damp 1 f_max 120 N_f 100 seed 35082 eta 1 beta 5 T_init 300
QBMSST parameters:
Shock in z direction
Cell mass-like parameter qmass (units of mass^2/length^4) = 4.00000e+01
Shock velocity = 7.80000e+01
Artificial viscosity (units of mass/length/time) = 0.00000e+00
Initial pressure calculated on first step
Initial volume calculated on first step
Initial energy calculated on first step
Resetting global fix info from restart file:
fix style: qbmsst, fix ID: shock
fix_modify shock energy yes
variable dhug equal f_shock[1]
variable dray equal f_shock[2]
variable lgr_vel equal f_shock[3]
variable lgr_pos equal f_shock[4]
variable T_qm equal f_shock[5] #Temperature with quantum nuclear correction
thermo_style custom step v_T_qm press etotal vol lx ly lz pzz v_dhug v_dray v_lgr_vel v_lgr_pos
thermo 500
timestep ${delta_t}
timestep 0.001
restart 1000 restart
run 10000 #10 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.302953
grid = 9 8 16
stencil order = 5
estimated absolute RMS force accuracy = 0.00105569
estimated relative force accuracy = 7.33134e-05
using double precision FFTW3
3d grid and FFT values/proc = 2640 288
All restart file global fix info was re-assigned
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 2 4
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
(2) pair table, perpetual, skip from (1)
attributes: half, newton on
pair build: skip
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 79.71 | 79.71 | 79.71 Mbytes
Step v_T_qm Press TotEng Volume Lx Ly Lz Pzz v_dhug v_dray v_lgr_vel v_lgr_pos
1000 281.98328 -6031.2395 -2778.6227 1881.8369 9.9267847 8.596848 22.051297 -14113.621 1.3373278 13060.059 -0.79055793 0
1500 266.12746 44405.573 -2777.9815 1739.6543 9.9267847 8.596848 20.385206 92590.239 -12.06041 397.47049 5.1624821 -37.823748
2000 255.79411 17620.408 -2777.9685 1785.7619 9.9267847 8.596848 20.925494 48670.364 -16.082827 -4813.6764 3.2320016 -73.974437
2500 256.8887 40153.833 -2778.4337 1752.9461 9.9267847 8.596848 20.540959 79665.002 7.7413878 -1368.8927 4.6059671 -112.35254
3000 261.55251 5315.4799 -2779.0755 1834.3375 9.9267847 8.596848 21.4947 15896.368 22.588205 3192.882 1.1981949 -148.36068
3500 261.57101 57911.809 -2778.1223 1713.3956 9.9267847 8.596848 20.077507 110996.8 -9.4471543 -3240.9018 6.2619064 -186.41261
4000 254.88665 13952.95 -2778.4816 1818.2782 9.9267847 8.596848 21.306518 26833.588 2.2818412 647.88057 1.8705799 -222.72504
4500 240.08908 73322.997 -2776.7382 1668.6666 9.9267847 8.596848 19.553375 151978.11 -43.917346 189.1572 8.1346613 -260.52885
5000 214.49084 1925.2557 -2777.0657 1890.0985 9.9267847 8.596848 22.148106 -5218.7292 -44.5537 28890.787 -1.1364617 -297.26329
5500 194.6515 71804.842 -2777.3417 1669.7297 9.9267847 8.596848 19.565832 146911.42 -34.911593 -3985.0635 8.0901523 -334.1879
6000 186.23814 10196.007 -2777.1394 1837.3793 9.9267847 8.596848 21.530344 23550.907 -18.381207 13401.096 1.0708382 -371.9208
6500 172.53603 5474.3725 -2777.4502 1818.0038 9.9267847 8.596848 21.303303 18389.825 -22.65951 -8026.2088 1.8820667 -407.83084
7000 160.91186 107908.64 -2777.6746 1621.7378 9.9267847 8.596848 19.003464 196841.27 -8.6606903 5654.1938 10.099523 -444.9925
7500 146.01905 147030.69 -2777.2543 1539.7536 9.9267847 8.596848 18.042777 253089.02 -43.928324 -6926.1018 13.532114 -478.63113
8000 207.17758 837859.1 -2796.8957 989.32874 9.9267847 8.596848 11.592918 811765.11 1172.3778 89652.363 36.577833 -503.41923
8500 725.15657 853732.89 -2832.3144 974.18299 9.9267847 8.596848 11.415441 773926.64 1749.5702 39098.598 37.21197 -524.17835
9000 1554.6089 807867.74 -2843.0063 990.10922 9.9267847 8.596848 11.602064 749697.22 1959.0322 28239.71 36.545155 -544.77354
9500 2440.1194 748145.05 -2839.2364 992.38871 9.9267847 8.596848 11.628775 691503.58 1437.0708 -28040.223 36.449715 -565.41198
10000 3112.1817 823862.43 -2820.0495 982.35471 9.9267847 8.596848 11.511197 754954.89 1330.6807 26987.244 36.869828 -586.12357
10500 3550.0273 868916.79 -2803.7678 983.70386 9.9267847 8.596848 11.527006 867368.45 1727.9058 140533.46 36.813341 -607.00946
11000 3839.7527 830581.55 -2795.3804 995.31485 9.9267847 8.596848 11.663063 811740 1150.0462 94652.768 36.327201 -628.02229
Loop time of 15.1476 on 4 procs for 10000 steps with 144 atoms
Performance: 57.039 ns/day, 0.421 hours/ns, 660.171 timesteps/s
91.3% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.7228 | 9.085 | 10.626 | 36.0 | 59.98
Kspace | 1.6343 | 3.1795 | 4.5467 | 61.0 | 20.99
Neigh | 0.02063 | 0.027076 | 0.034395 | 3.1 | 0.18
Comm | 0.54719 | 0.57781 | 0.60468 | 2.8 | 3.81
Output | 0.10128 | 0.1019 | 0.10373 | 0.3 | 0.67
Modify | 2.0819 | 2.1159 | 2.1495 | 1.8 | 13.97
Other | | 0.06035 | | | 0.40
Nlocal: 36 ave 38 max 33 min
Histogram: 1 0 0 0 0 0 1 0 1 1
Nghost: 4267 ave 4304 max 4239 min
Histogram: 1 0 1 0 1 0 0 0 0 1
Neighs: 18859.2 ave 25108 max 12333 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Total # of neighbors = 75437
Ave neighs/atom = 523.868
Neighbor list builds = 95
Dangerous builds = 0
Total wall time: 0:00:17

2
examples/USER/qtb/alpha_quartz_qbmsst/potential_SiO2.TPF Executable file → Normal file
View File

@ -1,4 +1,4 @@
#Potential between O,Si atoms. r ranges from 0.1 to 20.0. This file+BKS(with no shift)=alpha*(r-rs)^2(0.1<r<rs);BKS(20.0>r>rs). alpha=100, rsOO=1.43869, rsSiO=1.19362. All with metal units.
#Potential between O,Si atoms. UNITS: metal r ranges from 0.1 to 20.0. This file+BKS(with no shift)=alpha*(r-rs)^2(0.1<r<rs);BKS(20.0>r>rs). alpha=100, rsOO=1.43869, rsSiO=1.19362. All with metal units.
Si-O
N 39901

79
examples/USER/qtb/alpha_quartz_qtb/alpha_quartz_qtb.in
View File

@ -1,79 +0,0 @@
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 100
run 20000 # 20 ps
unfix quartz_qtb
unfix scapegoat_qtb

79
examples/USER/qtb/alpha_quartz_qtb/in.alpha_quartz_qtb Normal file
View File

@ -0,0 +1,79 @@
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 &
a1 4.916000 0.000000 0.000000 &
a2 -2.45800 4.257381 0.000000 &
a3 0.000000 0.000000 5.405400 &
&
basis 0.469700 0.000000 0.000000 &
basis 0.000000 0.469700 0.666667 &
basis 0.530300 0.530300 0.333333 &
&
basis 0.413500 0.266900 0.119100 &
basis 0.266900 0.413500 0.547567 &
basis 0.733100 0.146600 0.785767 &
basis 0.586500 0.853400 0.214233 &
basis 0.853400 0.586500 0.452433 &
basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
create_atoms 1 box &
basis 1 1 &
basis 2 1 &
basis 3 1 &
basis 4 2 &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2 &
basis 9 2
replicate ${x_rep} ${y_rep} ${z_rep}
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
set type 2 charge -1.2
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
unfix quartz_qtb
unfix scapegoat_qtb

152
examples/USER/qtb/alpha_quartz_qtb/log.15Jun20.alpha_quartz_qtb.g++.1 Normal file
View File

@ -0,0 +1,152 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 1 by 1 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 2 ${z_rep}
replicate 2 2 2
orthogonal box = (0.0 0.0 0.0) to (9.832 17.029524 10.8108)
1 by 1 by 1 MPI processor grid
144 atoms
replicate CPU = 0.000219584 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_style buck/coul/long 10
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.307414
grid = 9 15 10
stencil order = 5
estimated absolute RMS force accuracy = 0.000822922
estimated relative force accuracy = 5.71487e-05
using double precision FFTW3
3d grid and FFT values/proc = 5984 1350
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 3 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 79.54 | 79.54 | 79.54 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34025.794 -2793.6041 1810.0985 9.832 17.029524 10.8108 -37478.502 -37477.413 -27121.466 -1.3649088e-10 1.3388978 5.8209479e-10
500 281.29079 -40385.348 -2786.6755 1844.5575 9.893999 17.136909 10.878971 -44649.574 -45631.516 -30874.953 -5970.3691 3630.1324 5208.8966
1000 405.39848 -15491.657 -2783.3315 1874.0851 9.9465141 17.227868 10.936714 -7770.4561 -21469.887 -17234.627 -4706.5632 -8313.9522 -5109.7918
1500 428.98568 -12118.951 -2781.3232 1874.1627 9.9466513 17.228106 10.936865 -11239.135 -11740.052 -13377.666 -3778.9317 -6220.1431 12775.412
2000 522.11905 -6687.482 -2779.2181 1855.9626 9.914349 17.172157 10.901347 -8016.0133 -10737.23 -1309.2028 -4980.3805 5270.2674 5848.5479
2500 496.74376 4264.2623 -2778.9979 1849.3244 9.9025147 17.151659 10.888335 -477.1374 3487.19 9782.7343 -4961.2016 2380.6522 4736.0758
3000 456.49628 2320.781 -2779.3844 1853.2925 9.9095923 17.163918 10.896117 5479.6232 -2954.3023 4437.022 3719.9287 4445.0723 -3278.5058
3500 485.20722 -7480.1789 -2778.6062 1859.6305 9.920876 17.183462 10.908524 -9340.6334 -6129.8494 -6970.0541 -7379.3507 1772.8159 334.33057
4000 527.61216 -13499.73 -2777.3433 1889.9405 9.9744857 17.276316 10.96747 -16483.038 -7465.2297 -16550.923 -2517.02 -1863.063 3314.927
4500 519.94117 721.60614 -2777.8506 1879.6562 9.9563603 17.244922 10.947541 -913.2791 -1765.7541 4843.8516 4466.5704 -14141.087 -6439.5669
5000 505.27757 -6278.3805 -2777.3641 1881.2931 9.9592497 17.249927 10.950718 -14254.233 -2653.6233 -1927.2858 1838.1568 5767.9267 597.47761
5500 500.70903 11303.821 -2777.8881 1871.0076 9.9410666 17.218433 10.930724 -6452.7947 24876.967 15487.29 522.01171 10473.257 9780.893
6000 526.65329 7991.2419 -2777.172 1856.9227 9.9160583 17.175117 10.903227 -68.823156 11005.468 13037.081 1253.9214 10039.559 1053.0486
6500 485.30026 12811.546 -2777.5866 1845.31 9.8953442 17.139239 10.88045 10063.921 20215.037 8155.6798 -3886.954 2249.2807 4855.0011
7000 507.85472 2649.7919 -2777.3359 1861.2877 9.923822 17.188564 10.911763 -4214.7779 6995.1472 5169.0064 -2188.489 6157.0955 533.65478
7500 528.5729 3161.4629 -2779.0851 1855.7946 9.9140499 17.171639 10.901018 2935.365 -2873.1363 9422.1601 771.1885 -4360.9131 4939.8209
8000 533.77283 4534.849 -2777.6538 1858.4772 9.9188246 17.179909 10.906268 -1187.9433 15739.396 -946.90551 -5187.8588 2446.5059 8079.2032
8500 518.71765 1108.9877 -2777.7019 1866.6125 9.9332765 17.20494 10.922159 8720.4976 -8234.9325 2841.3979 5148.5004 -2125.3524 -4127.7468
9000 536.71495 -496.88283 -2778.0262 1877.7099 9.9529227 17.238968 10.943761 -3481.5874 -4611.6246 6602.5634 -2788.5111 -13323.148 4338.813
9500 527.06773 -236.09043 -2778.1125 1895.9227 9.9849986 17.294525 10.97903 -12233.409 7578.0514 3947.0863 -6399.0254 995.22838 8590.7109
10000 526.77335 -4480.6866 -2777.7171 1886.8998 9.9691335 17.267046 10.961585 -3139.961 1336.993 -11639.092 13496.371 -11543.676 -6180.9262
Loop time of 8.86837 on 1 procs for 10000 steps with 144 atoms
Performance: 97.425 ns/day, 0.246 hours/ns, 1127.603 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 6.1503 | 6.1503 | 6.1503 | 0.0 | 69.35
Kspace | 1.1522 | 1.1522 | 1.1522 | 0.0 | 12.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.11705 | 0.11705 | 0.11705 | 0.0 | 1.32
Output | 0.00035834 | 0.00035834 | 0.00035834 | 0.0 | 0.00
Modify | 1.4245 | 1.4245 | 1.4245 | 0.0 | 16.06
Other | | 0.02397 | | | 0.27
Nlocal: 144 ave 144 max 144 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 3804 ave 3804 max 3804 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 41952 ave 41952 max 41952 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
Total wall time: 0:00:08

152
examples/USER/qtb/alpha_quartz_qtb/log.15Jun20.alpha_quartz_qtb.g++.4 Normal file
View File

@ -0,0 +1,152 @@
LAMMPS (15 Jun 2020)
using 1 OpenMP thread(s) per MPI task
## This script first constructs an alpha quartz structure of a given size. It then uses fix qtb to equilibrate the computational cell to the specified temperature and pressure.
variable x_rep equal 2 #x-direction replication number
variable y_rep equal 2 #y-direction replication number
variable z_rep equal 2 #z-direction replication number
variable cut_off equal 10.0 #Cut-off distance for the Buckingham term (Angstrom in metal units)
variable pressure equal 1.03125 #Initial state pressure (bar in metal units)
variable temperature equal 300.0 #Initial state quantum temperature (K in metal units)
variable delta_t equal 1.0e-3 #MD timestep length (ps in metal units)
variable damp_qtb equal 1 #1/gamma where gamma is the friction coefficient in quantum thermal bath (ps in metal units)
## This part defines units, alpha-quartz crystal, and atomic information
#General
units metal
dimension 3
boundary p p p
atom_style charge
#Lattice
lattice custom 1.0 a1 4.916000 0.000000 0.000000 a2 -2.45800 4.257381 0.000000 a3 0.000000 0.000000 5.405400 basis 0.469700 0.000000 0.000000 basis 0.000000 0.469700 0.666667 basis 0.530300 0.530300 0.333333 basis 0.413500 0.266900 0.119100 basis 0.266900 0.413500 0.547567 basis 0.733100 0.146600 0.785767 basis 0.586500 0.853400 0.214233 basis 0.853400 0.586500 0.452433 basis 0.146600 0.733100 0.880900 #American Mineralogist 65 920 1980 (Space Group 154)
Lattice spacing in x,y,z = 7.374 4.25738 5.4054
#Computational Cell
region orthorhombic_unit_cell block 0 4.916000 0 8.514762 0 5.405400 units box
create_box 2 orthorhombic_unit_cell
Created orthogonal box = (0.0 0.0 0.0) to (4.916 8.514762 5.4054)
1 by 2 by 2 MPI processor grid
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2
Created 18 atoms
create_atoms CPU = 0.000 seconds
replicate ${x_rep} ${y_rep} ${z_rep}
replicate 2 ${y_rep} ${z_rep}
replicate 2 2 ${z_rep}
replicate 2 2 2
orthogonal box = (0.0 0.0 0.0) to (9.832 17.029524 10.8108)
1 by 2 by 2 MPI processor grid
144 atoms
replicate CPU = 0.000231981 secs
#Atomic Information
mass 1 28.085500
mass 2 15.999400
set type 1 charge +2.4
48 settings made for charge
set type 2 charge -1.2
96 settings made for charge
## This part implements the BKS pair potential with a cut-off distance for the Buckingham term. Long range Coulomb interactions are evaluated with the pppm method.
#Pair Potentials
pair_style buck/coul/long ${cut_off} #BKS interaction, PRL 64 1955 (1990)
pair_style buck/coul/long 10
pair_coeff 1 1 0.0 1.0 0.0
pair_coeff 1 2 18003.757200 0.205205 133.538100
pair_coeff 2 2 1388.773000 0.362319 175.000000
pair_modify shift yes
kspace_style pppm 1.0e-4
#Neighbor style
neighbor 2.0 bin
neigh_modify check yes every 1 delay 0 page 100000 one 2000
## This part equilibrates your crystal to a pressure of ${pressure}(unit pressure) and a temperature of ${temperature}(unit temperatureture) with quantum nuclear effects
variable p_damp equal ${delta_t}*1000 #Recommended pressure damping parameter in fix nph
variable p_damp equal 0.001*1000
fix scapegoat_qtb all nph iso ${pressure} ${pressure} ${p_damp} ptemp ${temperature} #NPH does the time integration
fix scapegoat_qtb all nph iso 1.03125 ${pressure} ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 ${p_damp} ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp ${temperature}
fix scapegoat_qtb all nph iso 1.03125 1.03125 1 ptemp 300
fix quartz_qtb all qtb temp ${temperature} damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100 #Change f_max (THz) if your Debye frequency is higher
fix quartz_qtb all qtb temp 300 damp ${damp_qtb} seed 35082 f_max 120.00 N_f 100
fix quartz_qtb all qtb temp 300 damp 1 seed 35082 f_max 120.00 N_f 100
thermo_style custom step temp press etotal vol lx ly lz pxx pyy pzz pxy pyz pxz
thermo 500
run 10000 # 20 ps
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
G vector (1/distance) = 0.307414
grid = 9 15 10
stencil order = 5
estimated absolute RMS force accuracy = 0.000822922
estimated relative force accuracy = 5.71487e-05
using double precision FFTW3
3d grid and FFT values/proc = 2688 405
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 2 3 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair buck/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 79.08 | 79.08 | 79.09 Mbytes
Step Temp Press TotEng Volume Lx Ly Lz Pxx Pyy Pzz Pxy Pyz Pxz
0 0 -34025.794 -2793.6041 1810.0985 9.832 17.029524 10.8108 -37478.502 -37477.413 -27121.466 -5.9958847e-10 1.3388978 7.2750373e-10
500 222.04947 -50221.579 -2787.6677 1851.5661 9.9065143 17.158586 10.892732 -61493.697 -53512.432 -35658.61 4973.9976 7095.5481 -2041.6341
1000 361.92367 -14345.85 -2783.1509 1861.0579 9.9234137 17.187857 10.911314 -4145.4149 -28701.195 -10190.939 7896.3934 -3901.2874 -490.57107
1500 457.97039 -4934.8727 -2779.8321 1860.2254 9.9219337 17.185294 10.909687 -3680.7192 -3045.0707 -8078.8283 456.70383 -4705.3346 -2175.8144
2000 523.52684 -9516.755 -2778.5181 1866.3577 9.9328244 17.204157 10.921662 -11042.489 -7777.5634 -9730.2127 -2016.3336 6027.001 -4150.3656
2500 489.58881 -4968.5157 -2777.3948 1864.0745 9.9287723 17.197139 10.917206 -13652.344 -2823.514 1570.3111 -7481.9537 -1150.3548 10502.368
3000 559.52782 -2882.7076 -2777.5527 1883.2223 9.9626528 17.255821 10.95446 3061.1755 -10570.656 -1138.642 -12045.354 -856.20951 16292.443
3500 521.67929 9974.5929 -2776.7752 1880.5936 9.9580152 17.247789 10.94936 15358.559 7855.8683 6709.3509 7292.9372 -9848.9204 -523.61056
4000 497.90872 -2012.9259 -2776.5554 1862.1703 9.9253904 17.191281 10.913488 -1154.5123 4270.0275 -9154.2927 971.94826 -10157.618 4694.0509
4500 533.64016 -7218.9278 -2775.8789 1883.3041 9.962797 17.256071 10.954618 -18299.547 -5497.566 2140.3296 -1335.6063 -10353.21 5703.7506
5000 551.61416 1590.9702 -2777.6093 1866.2047 9.9325531 17.203687 10.921363 -4600.02 6535.3 2837.6306 3412.3383 9492.18 1017.5742
5500 499.36075 188.82067 -2777.9872 1863.2925 9.9273838 17.194734 10.91568 -3238.914 1143.013 2662.363 4193.7623 -11565.423 2575.9361
6000 478.563 4064.8319 -2778.946 1867.7185 9.935238 17.208337 10.924316 1947.7246 3346.7411 6900.0301 -6339.9554 4133.6942 -4555.406
6500 512.63865 10227.461 -2778.8476 1855.5323 9.9135828 17.17083 10.900505 7423.8967 7558.2024 15700.285 -621.4585 -2620.4837 -3256.7524
7000 489.9889 13037.303 -2778.8793 1856.2469 9.9148553 17.173034 10.901904 10690.345 16770.786 11650.779 -4056.5527 -5023.8847 469.21909
7500 495.52187 1320.5068 -2778.1189 1871.7467 9.9423755 17.2207 10.932164 1978.2905 738.78041 1244.4496 1826.0923 -7829.3563 1873.2713
8000 474.60945 -4203.2068 -2778.8915 1866.5966 9.9332482 17.204891 10.922128 -1480.6896 -12516.261 1387.3306 2731.4462 -1292.9741 10743.939
8500 473.16225 -6266.1992 -2778.594 1872.9075 9.9444304 17.224259 10.934423 -12680.492 -2832.6603 -3285.4455 7226.9632 3762.6841 -5834.9064
9000 486.6579 2843.7947 -2778.0388 1877.3735 9.9523282 17.237939 10.943107 805.23659 6213.7247 1512.4228 2685.2063 -3517.5266 -17054.035
9500 549.35112 -1028.3899 -2776.8124 1880.7965 9.9583733 17.248409 10.949754 -1817.8413 2754.8459 -4022.1743 -3101.1463 8397.2345 -8608.1342
10000 562.27081 12885.53 -2775.7435 1850.2864 9.9042316 17.154633 10.890222 15758.218 9989.5121 12908.859 -25.724137 -16691.374 267.85371
Loop time of 3.80648 on 4 procs for 10000 steps with 144 atoms
Performance: 226.981 ns/day, 0.106 hours/ns, 2627.100 timesteps/s
94.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.3526 | 1.6581 | 1.9634 | 21.3 | 43.56
Kspace | 0.92143 | 1.2222 | 1.5232 | 24.7 | 32.11
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.31886 | 0.32256 | 0.32604 | 0.5 | 8.47
Output | 0.0003643 | 0.00083923 | 0.0022533 | 0.0 | 0.02
Modify | 0.39166 | 0.45985 | 0.52607 | 8.9 | 12.08
Other | | 0.143 | | | 3.76
Nlocal: 36 ave 36 max 36 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 2614 ave 2614 max 2614 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 10488 ave 12570 max 8406 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 41952
Ave neighs/atom = 291.333
Neighbor list builds = 0
Dangerous builds = 0
unfix quartz_qtb
unfix scapegoat_qtb
Total wall time: 0:00:03

18
examples/threebody/Si.sw
View File

@ -1,18 +0,0 @@
# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985)
# Stillinger-Weber parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units:
# epsilon = eV; sigma = Angstroms
# other quantities are unitless
# format of a single entry (one or more lines):
# element 1, element 2, element 3,
# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol
# Here are the original parameters in metal units, for Silicon from:
#
# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985)
#
Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333
7.049556277 0.6022245584 4.0 0.0 0.0

1
examples/threebody/Si.sw Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Si.sw

71
lib/awpmd/ivutils/include/pairhash.h
View File

@ -122,12 +122,9 @@
# define PAIRHASH_H
/*e @file pairhash.h @brief pair hash table
*/
*/
/*r @file pairhash.h @brief ðàáîòà ñ õåø-òàáëèöàìè ïàðíûõ âåëè÷èí
*/
# include "refobj.h"
@ -146,7 +143,7 @@ public:
incr=inc_first ? parent->sizex : 1 ;
};
iterator(T *ptr_,size_t incr_):ptr(ptr_),incr(incr_){}
public:
public:
iterator(const iterator &other):ptr(other.ptr),incr(other.incr){
}
iterator():ptr(NULL),incr(0){}
@ -159,7 +156,7 @@ public:
++*this;
return tmp;
}
iterator operator+(int delta) const {
iterator operator+(int delta) const {
return iterator(ptr+delta*incr,incr);
}
bool operator!=(const iterator &other) const {
@ -194,7 +191,7 @@ public:
init(other.sizex,other.sizey,1);
size_t n=get_datasize(sizex,sizey);
for(size_t i=0;i<n;i++)
arr[i]=other.arr[i];
arr[i]=other.arr[i];
}
return *this;
}
@ -224,7 +221,7 @@ public:
virtual int init(size_t nx, size_t ny, int smanaged=-1){
int managed=parr.managed();
if(managed && (sizex!=nx || sizey!=ny)){
if(managed && (sizex!=nx || sizey!=ny)){
parr.reset(NULL,0);
}
if(smanaged>=0){ // for changing the managed flag?
@ -241,14 +238,14 @@ public:
arr=parr.ptr();
}
return 1;
}
}
recmatrix(size_t nx, size_t ny):sizex(0), sizey(0){
recmatrix(size_t nx, size_t ny):sizex(0), sizey(0){
init(nx,ny,1);
}
//e initializes by unmanaged pointer
recmatrix(size_t nx, size_t ny , T *ptr):parr(ptr,0),sizex(nx), sizey(ny) {
recmatrix(size_t nx, size_t ny , T *ptr):parr(ptr,0),sizex(nx), sizey(ny) {
init(nx,ny);
}
@ -264,14 +261,14 @@ public:
size_t i, n=get_datasize(sizex,sizey);
for(i=0;i<n;i++)arr[i]=val;
}
void SetDiag(const T &val){
size_t i, size=(sizex>sizey? sizey: sizex);
for(i=0;i<size;i++){
(*this)(i,i)=val;
}
}
/// returns iterator with fixed first index to iterate through matrix line elements
iterator fix_first(size_t first, size_t second) const {
return iterator(this,first,second,false);
@ -281,14 +278,14 @@ public:
iterator fix_second(size_t first, size_t second) const {
return iterator(this,first,second, true);
}
};
//e square matrix
template <class T>
class sqmatrix: public recmatrix<T> {
public:
size_t size;
@ -313,18 +310,18 @@ public:
return n*n;
}
virtual int init(size_t n, int smanaged=-1){
size=n;
return recmatrix<T>::init(n,n,smanaged);
}
}
sqmatrix(size_t n):size(0){
sqmatrix(size_t n):size(0){
init(n,1);
}
//e initializes by unmanaged pointer
sqmatrix(size_t n, T * /* ptr */):size(n){
sqmatrix(size_t n, T * /* ptr */):size(n){
init(n);
}
@ -335,7 +332,7 @@ public:
recmatrix<T>::parr.reset(ptr,0);
}
};
@ -355,7 +352,7 @@ public:
incr=inc_first ? parent->size : 1 ;
};
iterator(T *ptr_,size_t incr_):ptr(ptr_),incr(incr_){}
public:
public:
iterator(const iterator &other):ptr(other.ptr),incr(other.incr){
}
iterator():ptr(NULL),incr(0){}
@ -368,7 +365,7 @@ public:
++*this;
return tmp;
}
iterator operator+(int delta) const {
iterator operator+(int delta) const {
return iterator(ptr+delta*incr,incr);
}
bool operator!=(const iterator &other) const {
@ -403,7 +400,7 @@ public:
init(other.size,1);
size_t n=get_datasize(size);
for(size_t i=0;i<n;i++)
arr[i]=other.arr[i];
arr[i]=other.arr[i];
}
return *this;
}
@ -432,7 +429,7 @@ public:
virtual int init(size_t n, int smanaged=-1){
int managed=parr.managed();
if(managed && size!=n){
if(managed && size!=n){
parr.reset(NULL,0);
}
if(smanaged>=0){ // for changing the managed flag?
@ -448,14 +445,14 @@ public:
arr=parr.ptr();
}
return 1;
}
}
sqmatrix(size_t n):size(0){
sqmatrix(size_t n):size(0){
init(n,1);
}
//e initializes by unmanaged pointer
sqmatrix(size_t n, T *ptr):parr(ptr,0),size(n){
sqmatrix(size_t n, T *ptr):parr(ptr,0),size(n){
init(n);
}
@ -470,14 +467,14 @@ public:
size_t i, n=get_datasize(size);
for(i=0;i<n;i++)arr[i]=val;
}
void SetDiag(const T &val){
size_t i;
for(i=0;i<size;i++){
(*this)(i,i)=val;
}
}
/// returns iterator with fixed first index to iterate through matrix line elements
iterator fix_first(size_t first, size_t second) const {
return iterator(this,first,second,false);
@ -487,15 +484,15 @@ public:
iterator fix_second(size_t first, size_t second) const {
return iterator(this,first,second, true);
}
};
# endif
# endif
//e prints the matrix into a file
template< class matrix_t>
int fileout(FILE *f, const matrix_t &matr, const char *elm_fmt, const char *elm_sep=" ", const char *line_sep="\n"){
size_t i, j;
size_t i, j;
int res=0;
for(i=0;i<matr.size;i++){
for(j=0;j<matr.size;j++){
@ -513,15 +510,15 @@ template <class T>
class smatrix: public sqmatrix<T>{
typedef sqmatrix<T> base_t;
public:
virtual size_t get_datasize(size_t n) const{
return n*(n+1)/2;
}
size_t index(size_t i, size_t j) const {
if(i>=j)
return (2*base_t::size-j-1)*j/2+i;
else
else
return (2*base_t::size-i-1)*i/2+j;
}
@ -572,7 +569,7 @@ public:
hmatrix() : smatrix<T>() {}
hmatrix(const smatrix<T> &other) : smatrix<T>(other) {}
//e copy constructor: makes a managed copy
hmatrix(const hmatrix &other): smatrix<T>(other){}
@ -671,7 +668,7 @@ public:
indm.Set(-1);
return 1;
}
virtual ~PairHashM(){
delete [] arr;
}

16
lib/gpu/Makefile.hip
View File

@ -33,10 +33,10 @@ HIP_HOST_OPTS += -DMPI_GERYON -DUCL_NO_EXIT
# this settings should match LAMMPS Makefile
MPI_COMP_OPTS = $(shell mpicxx --showme:compile)
MPI_LINK_OPTS = $(shell mpicxx --showme:link)
#MPI_COMP_OPTS += -I/usr/include/mpi -DMPICH_IGNORE_CXX_SEEK -DOMPI_SKIP_MPICXX=1
HIP_PATH ?= $(wildcard /opt/rocm/hip)
HIP_PLATFORM=$(shell $(HIP_PATH)/bin/hipconfig --compiler)
HIP_PLATFORM=$(shell $(HIP_PATH)/bin/hipconfig --platform)
HIP_COMPILER=$(shell $(HIP_PATH)/bin/hipconfig --compiler)
ifeq (hcc,$(HIP_PLATFORM))
HIP_OPTS += -ffast-math
@ -48,8 +48,6 @@ else ifeq (nvcc,$(HIP_PLATFORM))
-gencode arch=compute_50,code=[sm_50,compute_50] -gencode arch=compute_52,code=[sm_52,compute_52] -gencode arch=compute_53,code=[sm_53,compute_53]\
-gencode arch=compute_60,code=[sm_60,compute_60] -gencode arch=compute_61,code=[sm_61,compute_61] -gencode arch=compute_62,code=[sm_62,compute_62]\
-gencode arch=compute_70,code=[sm_70,compute_70] -gencode arch=compute_72,code=[sm_72,compute_72] -gencode arch=compute_75,code=[sm_75,compute_75]
else
$(error Specify HIP platform using 'export HIP_PLATFORM=(hcc,nvcc)')
endif
BIN_DIR = .
@ -66,7 +64,15 @@ BSH = /bin/sh
HIP_OPTS += -DUSE_HIP $(HIP_PRECISION)
HIP_GPU_OPTS += $(HIP_OPTS) -I./
ifeq (hcc,$(HIP_PLATFORM))
ifeq (clang,$(HIP_COMPILER))
HIP_HOST_OPTS += -fPIC
HIP_GPU_CC = $(HIP_PATH)/bin/hipcc --genco
HIP_GPU_OPTS_S = --offload-arch=$(HIP_ARCH)
HIP_GPU_OPTS_E =
HIP_KERNEL_SUFFIX = .cpp
HIP_LIBS_TARGET = export HCC_AMDGPU_TARGET := $(HIP_ARCH)
export HCC_AMDGPU_TARGET := $(HIP_ARCH)
else ifeq (hcc,$(HIP_COMPILER))
HIP_HOST_OPTS += -fPIC
HIP_GPU_CC = $(HIP_PATH)/bin/hipcc --genco
HIP_GPU_OPTS_S = -t="$(HIP_ARCH)" -f=\"

18
lib/gpu/lal_lj_tip4p_long.cpp
View File

@ -250,22 +250,22 @@ void LJTIP4PLongT::copy_relations_data(int n, tagint *tag, int *map_array,
if (ago == 0) {
hneight.zero();
{
UCL_H_Vec<tagint> host_tag_write(nall,*(this->ucl_device),UCL_WRITE_ONLY);
UCL_H_Vec<tagint> host_tag_write;
host_tag_write.view(tag, nall, *(this->ucl_device));
this->tag.resize_ib(nall);
for(int i=0; i<nall; ++i) host_tag_write[i] = tag[i];
ucl_copy(this->tag, host_tag_write, nall, false);
ucl_copy(this->tag, host_tag_write, false);
}
UCL_H_Vec<int> host_write(max_same,*(this->ucl_device),UCL_WRITE_ONLY);
UCL_H_Vec<int> host_write;
host_write.view(sametag, max_same, *(this->ucl_device));
this->atom_sametag.resize_ib(max_same);
for(int i=0; i<max_same; ++i) host_write[i] = sametag[i];
ucl_copy(this->atom_sametag, host_write, max_same, false);
ucl_copy(this->atom_sametag, host_write, false);
host_write.resize_ib(map_size);
this->map_array.resize_ib(map_size);
for(int i=0; i<map_size; ++i) host_write[i] = map_array[i];
ucl_copy(this->map_array, host_write, map_size, false);
host_write.view(map_array, map_size, *(this->ucl_device));
ucl_copy(this->map_array, host_write, false);
}
}

46
lib/gpu/lal_lj_tip4p_long.cu
View File

@ -110,10 +110,8 @@ __kernel void k_lj_tip4p_long_distrib(const __global numtyp4 *restrict x_,
const int typeO, const int typeH,
const numtyp alpha,
const __global numtyp *restrict q_, const __global acctyp4 *restrict ansO) {
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+tid;
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+THREAD_ID_X;
acctyp4 f;
f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
@ -122,6 +120,8 @@ __kernel void k_lj_tip4p_long_distrib(const __global numtyp4 *restrict x_,
int itype = ix.w;
acctyp4 fM, vM;
acctyp eM;
// placement of the virial in engv depends on eflag value
int engv_iter = eflag ? 2 : 0;
if (itype == typeH) {
int iO = hneigh[i*4];
if (iO < inum) {
@ -131,13 +131,13 @@ __kernel void k_lj_tip4p_long_distrib(const __global numtyp4 *restrict x_,
f.z += fM.z * (acctyp)0.5 * alpha;
if (vflag > 0) {
vM = ansO[inum +iO];
engv[inum*2 + i] += vM.x * (acctyp)0.5 * alpha;
engv[inum*3 + i] += vM.y * (acctyp)0.5 * alpha;
engv[inum*4 + i] += vM.z * (acctyp)0.5 * alpha;
engv[inum*engv_iter + i] += vM.x * (acctyp)0.5 * alpha; engv_iter++;
engv[inum*engv_iter + i] += vM.y * (acctyp)0.5 * alpha; engv_iter++;
engv[inum*engv_iter + i] += vM.z * (acctyp)0.5 * alpha; engv_iter++;
vM = ansO[inum*2+iO];
engv[inum*5 + i] += vM.x * (acctyp)0.5 * alpha;
engv[inum*6 + i] += vM.y * (acctyp)0.5 * alpha;
engv[inum*7 + i] += vM.z * (acctyp)0.5 * alpha;
engv[inum*engv_iter + i] += vM.x * (acctyp)0.5 * alpha; engv_iter++;
engv[inum*engv_iter + i] += vM.y * (acctyp)0.5 * alpha; engv_iter++;
engv[inum*engv_iter + i] += vM.z * (acctyp)0.5 * alpha;
}
}
} else {
@ -155,13 +155,13 @@ __kernel void k_lj_tip4p_long_distrib(const __global numtyp4 *restrict x_,
}
if (vflag > 0) {
vM = ansO[inum + i];
engv[inum*2 + i] += vM.x * (acctyp)(1 - alpha);
engv[inum*3 + i] += vM.y * (acctyp)(1 - alpha);
engv[inum*4 + i] += vM.z * (acctyp)(1 - alpha);
engv[inum*engv_iter + i] += vM.x * (acctyp)(1 - alpha); engv_iter++;
engv[inum*engv_iter + i] += vM.y * (acctyp)(1 - alpha); engv_iter++;
engv[inum*engv_iter + i] += vM.z * (acctyp)(1 - alpha); engv_iter++;
vM = ansO[inum*2 + i];
engv[inum*5 + i] += vM.x * (acctyp)(1 - alpha);
engv[inum*6 + i] += vM.y * (acctyp)(1 - alpha);
engv[inum*7 + i] += vM.z * (acctyp)(1 - alpha);
engv[inum*engv_iter + i] += vM.x * (acctyp)(1 - alpha); engv_iter++;
engv[inum*engv_iter + i] += vM.y * (acctyp)(1 - alpha); engv_iter++;
engv[inum*engv_iter + i] += vM.z * (acctyp)(1 - alpha);
}
}
acctyp4 old=ans[i];
@ -182,9 +182,8 @@ __kernel void k_lj_tip4p_reneigh(const __global numtyp4 *restrict x_,
const int typeO, const int typeH,
const __global tagint *restrict tag, const __global int *restrict map,
const __global int *restrict sametag) {
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+tid;
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+THREAD_ID_X;
if (i<nall) {
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
@ -237,13 +236,10 @@ __kernel void k_lj_tip4p_newsite(const __global numtyp4 *restrict x_,
__global numtyp4 *restrict m,
const int typeO, const int typeH,
const numtyp alpha, const __global numtyp *restrict q_) {
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+tid;
int i = BLOCK_ID_X*(BLOCK_SIZE_X)+THREAD_ID_X;
if (i<nall) {
int iO, iH1, iH2;
iO = i;
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
int itype = ix.w;
if (itype == typeO) {
@ -280,8 +276,6 @@ __kernel void k_lj_tip4p_long(const __global numtyp4 *restrict x_,
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
const numtyp eq_zero = 1e-6;
acctyp energy = (acctyp)0;
acctyp e_coul = (acctyp)0;
acctyp4 f, fO;
@ -595,8 +589,6 @@ __kernel void k_lj_tip4p_long_fast(const __global numtyp4 *restrict x_,
int tid, ii, offset;
atom_info(t_per_atom,ii,tid,offset);
const numtyp eq_zero = 1e-6;
__local numtyp4 lj1[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp4 lj3[MAX_SHARED_TYPES*MAX_SHARED_TYPES];
__local numtyp sp_lj[8];

2
potentials/MOH.nb3b.harmonic
View File

@ -1,4 +1,4 @@
# DATE: 2016-07-29 UNITS: metal CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none
# DATE: 2016-07-29 UNITS: real CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none
# nb3b/harmonic (nonbonded 3-body harmonic) parameters for various elements
#
# multiple entries can be added to this file, LAMMPS reads the ones it needs

2
potentials/charmm22.cmap
View File

@ -1,4 +1,4 @@
# DATE: 2016-09-26 CONTRIBUTOR: Robert Latour, latourr@clemson.edu CITATION: TBA
# DATE: 2016-09-26 UNITS: real CONTRIBUTOR: Robert Latour, latourr@clemson.edu CITATION: TBA
# Title: charmm22/charmm27 dihedral correction map
# alanine map, type 1

2
potentials/charmm36.cmap
View File

@ -1,4 +1,4 @@
# DATE: 2016-09-26 CONTRIBUTOR: Robert Latour, latourr@clemson.edu CITATION: TBA
# DATE: 2016-09-26 UNITS:real CONTRIBUTOR: Robert Latour, latourr@clemson.edu CITATION: TBA
# Title: charmm36 dihedral correction map
# alanine map, type 1

31
src/ASPHERE/fix_nph_asphere.cpp
View File

@ -13,6 +13,7 @@
#include "fix_nph_asphere.h"
#include <cstring>
#include <string>
#include "modify.h"
#include "error.h"
@ -35,35 +36,21 @@ FixNPHAsphere::FixNPHAsphere(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/asphere";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(tcmd + " all temp/asphere");
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
pcomputeflag = 1;
}

31
src/ASPHERE/fix_npt_asphere.cpp
View File

@ -13,6 +13,7 @@
#include "fix_npt_asphere.h"
#include <cstring>
#include <string>
#include "modify.h"
#include "error.h"
@ -34,35 +35,21 @@ FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/asphere";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(tcmd + " all temp/asphere");
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
pcomputeflag = 1;
}

18
src/ASPHERE/fix_nvt_asphere.cpp
View File

@ -13,9 +13,11 @@
#include "fix_nvt_asphere.h"
#include <cstring>
#include <string>
#include "group.h"
#include "modify.h"
#include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -33,17 +35,11 @@ FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string cmd = id + std::string("_temp");
id_temp = new char[cmd.size()+1];
strcpy(id_temp,cmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/asphere";
modify->add_compute(3,newarg);
delete [] newarg;
cmd += fmt::format(" {} temp/asphere",group->names[igroup]);
modify->add_compute(cmd);
tcomputeflag = 1;
}

31
src/BODY/fix_nph_body.cpp
View File

@ -17,6 +17,7 @@
#include "fix_nph_body.h"
#include <cstring>
#include <string>
#include "modify.h"
#include "error.h"
@ -38,35 +39,21 @@ FixNPHBody::FixNPHBody(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/body";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(tcmd + " all temp/body");
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
pcomputeflag = 1;
}

31
src/BODY/fix_npt_body.cpp
View File

@ -17,6 +17,7 @@
#include "fix_npt_body.h"
#include <cstring>
#include <string>
#include "modify.h"
#include "error.h"
@ -38,35 +39,21 @@ FixNPTBody::FixNPTBody(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/body";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(tcmd + " all temp/body");
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
pcomputeflag = 1;
}

16
src/BODY/fix_nvt_body.cpp
View File

@ -17,6 +17,7 @@
#include "fix_nvt_body.h"
#include <cstring>
#include <string>
#include "group.h"
#include "modify.h"
#include "error.h"
@ -37,17 +38,10 @@ FixNVTBody::FixNVTBody(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/body";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(tcmd + " all temp/body");
tcomputeflag = 1;
}

32
src/BODY/pair_body_rounded_polygon.cpp
View File

@ -236,7 +236,7 @@ void PairBodyRoundedPolygon::compute(int eflag, int vflag)
edge[jefirst+nj][4] = 0;
}
int interact, num_contacts, done;
int num_contacts, done;
double delta_a, j_a;
Contact contact_list[MAX_CONTACTS];
@ -244,15 +244,13 @@ void PairBodyRoundedPolygon::compute(int eflag, int vflag)
// check interaction between i's vertices and j' edges
interact = vertex_against_edge(i, j, k_nij, k_naij,
x, f, torque, tag, contact_list,
num_contacts, evdwl, facc);
vertex_against_edge(i, j, k_nij, k_naij, x, f, torque, tag,
contact_list, num_contacts, evdwl, facc);
// check interaction between j's vertices and i' edges
interact = vertex_against_edge(j, i, k_nij, k_naij,
x, f, torque, tag, contact_list,
num_contacts, evdwl, facc);
vertex_against_edge(j, i, k_nij, k_naij, x, f, torque, tag,
contact_list, num_contacts, evdwl, facc);
if (num_contacts >= 2) {
@ -595,7 +593,7 @@ void PairBodyRoundedPolygon::sphere_against_sphere(int i, int j,
{
double rradi,rradj;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double rij,rsqinv,R,fx,fy,fz,fpair,shift,energy;
double rij,rsqinv,R,fx,fy,fz,fn[3],ft[3],fpair,shift,energy;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
@ -641,6 +639,23 @@ void PairBodyRoundedPolygon::sphere_against_sphere(int i, int j,
vt1 = vr1 - vn1;
vt2 = vr2 - vn2;
vt3 = vr3 - vn3;
// normal friction term at contact
fn[0] = -c_n * vn1;
fn[1] = -c_n * vn2;
fn[2] = -c_n * vn3;
// tangential friction term at contact,
// excluding the tangential deformation term for now
ft[0] = -c_t * vt1;
ft[1] = -c_t * vt2;
ft[2] = -c_t * vt3;
fx += fn[0] + ft[0];
fy += fn[1] + ft[1];
fz += fn[2] + ft[2];
}
f[i][0] += fx;
@ -1349,4 +1364,3 @@ void PairBodyRoundedPolygon::distance(const double* x2, const double* x1,
+ (x2[1] - x1[1]) * (x2[1] - x1[1])
+ (x2[2] - x1[2]) * (x2[2] - x1[2]));
}

27
src/CLASS2/dihedral_class2.cpp
View File

@ -28,6 +28,7 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -207,24 +208,16 @@ void DihedralClass2::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
int me = comm->me;
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
error->warning(FLERR,fmt::format("Dihedral problem: {} {} {} {} {} {}",
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],
atom->tag[i3],atom->tag[i4]));
fmt::print(screen," 1st atom: {} {} {} {}\n",me,x[i1][0],x[i1][1],x[i1][2]);
fmt::print(screen," 2nd atom: {} {} {} {}\n",me,x[i2][0],x[i2][1],x[i2][2]);
fmt::print(screen," 3rd atom: {} {} {} {}\n",me,x[i3][0],x[i3][1],x[i3][2]);
fmt::print(screen," 4th atom: {} {} {} {}\n",me,x[i4][0],x[i4][1],x[i4][2]);
}
}

27
src/CLASS2/improper_class2.cpp
View File

@ -28,6 +28,7 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -153,24 +154,16 @@ void ImproperClass2::compute(int eflag, int vflag)
for (i = 0; i < 3; i++) {
if (costheta[i] == -1.0) {
int me;
MPI_Comm_rank(world,&me);
int me = comm->me;
if (screen) {
char str[128];
sprintf(str,"Improper problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
error->warning(FLERR,fmt::format("Improper problem: {} {} {} {} {} {}",
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],
atom->tag[i3],atom->tag[i4]));
fmt::print(screen," 1st atom: {} {} {} {}\n",me,x[i1][0],x[i1][1],x[i1][2]);
fmt::print(screen," 2nd atom: {} {} {} {}\n",me,x[i2][0],x[i2][1],x[i2][2]);
fmt::print(screen," 3rd atom: {} {} {} {}\n",me,x[i3][0],x[i3][1],x[i3][2]);
fmt::print(screen," 4th atom: {} {} {} {}\n",me,x[i4][0],x[i4][1],x[i4][2]);
}
}
}

20
src/CORESHELL/compute_temp_cs.cpp
View File

@ -19,6 +19,7 @@
#include "compute_temp_cs.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
@ -30,6 +31,7 @@
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -66,21 +68,13 @@ ComputeTempCS::ComputeTempCS(LAMMPS *lmp, int narg, char **arg) :
// create a new fix STORE style
// id = compute-ID + COMPUTE_STORE, fix group = compute group
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
id_fix = new char[n];
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
std::string fixcmd = id + std::string("_COMPUTE_STORE");
id_fix = new char[fixcmd.size()+1];
strcpy(id_fix,fixcmd.c_str());
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "0";
newarg[5] = (char *) "1";
modify->add_fix(6,newarg);
fixcmd += fmt::format(" {} STORE peratom 0 1", group->names[igroup]);
modify->add_fix(fixcmd);
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
// set fix store values = 0 for now
// fill them in via setup() once Comm::borders() has been called

16
src/GPU/pair_eam_alloy_gpu.cpp
View File

@ -371,8 +371,14 @@ void PairEAMAlloyGPU::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(PairEAM::lmp, filename,
"EAMAlloy", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -417,11 +423,19 @@ void PairEAMAlloyGPU::read_file(char *filename)
reader.next_dvector(&file->frho[i][1], file->nrho);
reader.next_dvector(&file->rhor[i][1], file->nr);
if (unit_convert) {
for (int j = 1; j < file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
}
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

16
src/GPU/pair_eam_fs_gpu.cpp
View File

@ -371,8 +371,14 @@ void PairEAMFSGPU::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(PairEAM::lmp, filename, "EAMFS",
unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -416,6 +422,10 @@ void PairEAMFSGPU::read_file(char *filename)
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
if (unit_convert) {
for (int j = 1; j <= file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
for (int j = 0; j < file->nelements; j++) {
reader.next_dvector(&file->rhor[i][j][1], file->nr);
@ -425,6 +435,10 @@ void PairEAMFSGPU::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k <= file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

5
src/GPU/pair_lj_cut_tip4p_long_gpu.cpp
View File

@ -109,10 +109,7 @@ PairLJCutTIP4PLongGPU::~PairLJCutTIP4PLongGPU()
void PairLJCutTIP4PLongGPU::compute(int eflag, int vflag)
{
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
ev_init(eflag,vflag);
int nall = atom->nlocal + atom->nghost;
int inum, host_start;

9
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -19,6 +19,7 @@
#include <mpi.h>
#include <cmath>
#include <cstring>
#include <string>
#include "atom.h"
#include "force.h"
#include "update.h"
@ -65,13 +66,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
// this is so final order of Modify:fix will conform to input script
fix_history = NULL;
char **fixarg = new char*[3];
fixarg[0] = (char *) "NEIGH_HISTORY_HH_DUMMY";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "DUMMY";
modify->add_fix(3,fixarg,1);
delete [] fixarg;
modify->add_fix("NEIGH_HISTORY_HH_DUMMY all DUMMY");
fix_dummy = (FixDummy *) modify->fix[modify->nfix-1];
}

9
src/GRANULAR/pair_granular.cpp
View File

@ -21,6 +21,7 @@
#include <mpi.h>
#include <cmath>
#include <cstring>
#include <string>
#include "atom.h"
#include "force.h"
#include "update.h"
@ -96,13 +97,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
// this is so final order of Modify:fix will conform to input script
fix_history = NULL;
char **fixarg = new char*[3];
fixarg[0] = (char *) "NEIGH_HISTORY_GRANULAR_DUMMY";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "DUMMY";
modify->add_fix(3,fixarg,1);
delete [] fixarg;
modify->add_fix("NEIGH_HISTORY_GRANULAR_DUMMY all DUMMY");
fix_dummy = (FixDummy *) modify->fix[modify->nfix-1];
}

6
src/KIM/kim_init.cpp
View File

@ -291,11 +291,7 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM
int ifix = modify->find_fix("KIM_MODEL_STORE");
if (ifix >= 0) modify->delete_fix(ifix);
char *fixarg[3];
fixarg[0] = (char *)"KIM_MODEL_STORE";
fixarg[1] = (char *)"all";
fixarg[2] = (char *)"STORE/KIM";
modify->add_fix(3,fixarg);
modify->add_fix("KIM_MODEL_STORE all STORE/KIM");
ifix = modify->find_fix("KIM_MODEL_STORE");
FixStoreKIM *fix_store = (FixStoreKIM *) modify->fix[ifix];

12
src/KOKKOS/atom_vec_kokkos.h
View File

@ -103,14 +103,14 @@ class AtomVecKokkos : public AtomVec {
ExecutionSpace space) = 0;
virtual int
pack_border_vel_kokkos(int n, DAT::tdual_int_2d k_sendlist,
DAT::tdual_xfloat_2d buf,int iswap,
int pbc_flag, int *pbc, ExecutionSpace space) { return 0; }
pack_border_vel_kokkos(int /*n*/, DAT::tdual_int_2d /*k_sendlist*/,
DAT::tdual_xfloat_2d /*buf*/,int /*iswap*/,
int /*pbc_flag*/, int * /*pbc*/, ExecutionSpace /*space*/) { return 0; }
virtual void
unpack_border_vel_kokkos(const int &n, const int &nfirst,
const DAT::tdual_xfloat_2d &buf,
ExecutionSpace space) {}
unpack_border_vel_kokkos(const int &/*n*/, const int & /*nfirst*/,
const DAT::tdual_xfloat_2d & /*buf*/,
ExecutionSpace /*space*/) {}
virtual int
pack_exchange_kokkos(const int &nsend, DAT::tdual_xfloat_2d &buf,

4
src/KOKKOS/kokkos_type.h
View File

@ -272,7 +272,7 @@ template<typename T1, typename T2>
class ScatterViewHelper<Kokkos::Experimental::ScatterDuplicated,T1,T2> {
public:
KOKKOS_INLINE_FUNCTION
static T1 get(const T1 &dup, const T2 &nondup) {
static T1 get(const T1 &dup, const T2 & /*nondup*/) {
return dup;
}
};
@ -1025,7 +1025,7 @@ struct params_lj_coul {
KOKKOS_INLINE_FUNCTION
params_lj_coul(){cut_ljsq=0;cut_coulsq=0;lj1=0;lj2=0;lj3=0;lj4=0;offset=0;};
KOKKOS_INLINE_FUNCTION
params_lj_coul(int i){cut_ljsq=0;cut_coulsq=0;lj1=0;lj2=0;lj3=0;lj4=0;offset=0;};
params_lj_coul(int /*i*/){cut_ljsq=0;cut_coulsq=0;lj1=0;lj2=0;lj3=0;lj4=0;offset=0;};
F_FLOAT cut_ljsq,cut_coulsq,lj1,lj2,lj3,lj4,offset;
};

15
src/KOKKOS/pair_eam_alloy_kokkos.cpp
View File

@ -985,8 +985,13 @@ void PairEAMAlloyKokkos<DeviceType>::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(lmp, filename, "EAMAlloy", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -1031,11 +1036,19 @@ void PairEAMAlloyKokkos<DeviceType>::read_file(char *filename)
reader.next_dvector(&file->frho[i][1], file->nrho);
reader.next_dvector(&file->rhor[i][1], file->nr);
if (unit_convert) {
for (int j = 1; j < file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
}
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

15
src/KOKKOS/pair_eam_fs_kokkos.cpp
View File

@ -985,8 +985,13 @@ void PairEAMFSKokkos<DeviceType>::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(lmp, filename, "EAMFS", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -1030,6 +1035,10 @@ void PairEAMFSKokkos<DeviceType>::read_file(char *filename)
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
if (unit_convert) {
for (int j = 1; j <= file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
for (int j = 0; j < file->nelements; j++) {
reader.next_dvector(&file->rhor[i][j][1], file->nr);
@ -1039,6 +1048,10 @@ void PairEAMFSKokkos<DeviceType>::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

4
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -632,8 +632,8 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeNeigh,const typen
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPBeta,const int& ii) const {
if (ii >= chunk_size) return;
if (ii >= chunk_size) return;
const int iatom_mod = ii % 32;
const int iatom_div = ii / 32;

6
src/KOKKOS/sna_kokkos.h
View File

@ -160,17 +160,17 @@ inline
t_sna_3c_ll ulist;
t_sna_3c_ll ylist;
// derivatives of data
t_sna_4c3_ll dulist;
// Modified structures for GPU backend
t_sna_3d_ll ulisttot_re; // split real,
t_sna_3d_ll ulisttot_im; // imag
t_sna_4c_ll ulisttot_pack; // AoSoA layout
t_sna_4c_ll zlist_pack; // AoSoA layout
t_sna_4d_ll blist_pack;
t_sna_4d_ll ylist_pack_re; // split real,
t_sna_4d_ll ylist_pack_re; // split real,
t_sna_4d_ll ylist_pack_im; // imag AoSoA layout
int idxcg_max, idxu_max, idxz_max, idxb_max;

37
src/KOKKOS/sna_kokkos_impl.h
View File

@ -295,7 +295,7 @@ void SNAKokkos<DeviceType>::pre_ui(const typename Kokkos::TeamPolicy<DeviceType>
// if m is on the "diagonal", initialize it with the self energy.
// Otherwise zero it out
double re_part = 0.;
double re_part = 0.;
if (m % (j+2) == 0 && (!chem_flag || ielem == jelem || wselfall_flag)) { re_part = wself; }
ulisttot_re(jjup, jelem, iatom) = re_part;
@ -368,7 +368,6 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
for (int j = 1; j <= twojmax; j++) {
const int jju = idxu_block[j];
const int jjup = idxu_block[j-1];
// fill in left side of matrix layer from previous layer
@ -436,7 +435,7 @@ void SNAKokkos<DeviceType>::compute_ui(const typename Kokkos::TeamPolicy<DeviceT
} else {
u_accum.re = -u_accum.re;
}
buf2[jju_shared_flip] = u_accum;
// split re, im to get fully coalesced atomic add
@ -552,10 +551,7 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
int idouble = 0;
for (int elem1 = 0; elem1 < nelements; elem1++) {
for (int elem2 = 0; elem2 < nelements; elem2++) {
const int jalloy = idouble;
for (int elem3 = 0; elem3 < nelements; elem3++) {
const int jjballoy = itriple;
double sumzu = 0.0;
double sumzu_temp = 0.0;
@ -566,7 +562,7 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
const int jjz_index = jjz+mb*(j+1)+ma;
if (2*mb == j) return; // I think we can remove this?
const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div);
sumzu_temp += utot.re * zloc.re + utot.im * zloc.im;
}
}
@ -582,7 +578,7 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma;
const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div);
sumzu_temp += utot.re * zloc.re + utot.im * zloc.im;
}
@ -593,7 +589,7 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma;
const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div);
const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div);
sumzu += 0.5 * (utot.re * zloc.re + utot.im * zloc.im);
} // end if jeven
@ -607,10 +603,10 @@ void SNAKokkos<DeviceType>::compute_bi(const int& iatom_mod, const int& jjb, con
sumzu -= bzero[j];
}
}
blist_pack(iatom_mod, jjb, jjballoy, iatom_div) = sumzu;
blist_pack(iatom_mod, jjb, itriple, iatom_div) = sumzu;
//} // end loop over j
//} // end loop over j1, j2
itriple++;
itriple++;
} // end loop over elem3
idouble++;
} // end loop over elem2
@ -791,7 +787,7 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
// single has a warp barrier at the end
Kokkos::single(Kokkos::PerThread(team), [=]() {
ulist_buf1[0] = {1., 0.};
dulist_buf1[0] = {0., 0.};
});
@ -881,7 +877,7 @@ void SNAKokkos<DeviceType>::compute_fused_deidrj(const typename Kokkos::TeamPoli
// u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])
if (j%2==1 && mb+1==n_mb) {
int sign_factor = (((ma+mb)%2==0)?1:-1);
const int jju_shared_flip = (j+1-mb)*(j+1)-(ma+1);
if (sign_factor == 1) {
@ -1078,12 +1074,11 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
int idouble = 0;
for (int elem1 = 0; elem1 < nelements; elem1++) {
for (int elem2 = 0; elem2 < nelements; elem2++) {
const auto jalloy = idouble;
auto jalloy = idouble; // must be non-const to work around gcc compiler bug
for (int elem3 = 0; elem3 < nelements; elem3++) {
Kokkos::parallel_for(Kokkos::TeamThreadRange(team,idxb_max),
[&] (const int& jjb) {
//for(int jjb = 0; jjb < idxb_max; jjb++) {
const auto jjballoy = itriple;
const int j1 = idxb(jjb, 0);
const int j2 = idxb(jjb, 1);
const int j = idxb(jjb, 2);
@ -1141,12 +1136,12 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
if (elem1 == elem2 && elem1 == elem3) {
sumzu -= bzero[j];
}
} else {
} else {
sumzu -= bzero[j];
}
}
blist(jjb, jjballoy, iatom) = sumzu;
blist(jjb, itriple, iatom) = sumzu;
});
});
//} // end loop over j
@ -1156,7 +1151,7 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
idouble++;
} // end loop over elem2
} // end loop over elem1
}
/* ----------------------------------------------------------------------
@ -2027,7 +2022,6 @@ double SNAKokkos<DeviceType>::memory_usage()
{
int jdimpq = twojmax + 2;
int jdim = twojmax + 1;
int natom_pad = ((natom + 32 - 1) / 32) * 32; // for AoSoA layouts
double bytes;
bytes = 0;
@ -2035,11 +2029,10 @@ double SNAKokkos<DeviceType>::memory_usage()
bytes += jdimpq*jdimpq * sizeof(double); // pqarray
bytes += idxcg_max * sizeof(double); // cglist
#ifdef KOKKOS_ENABLE_CUDA
if (std::is_same<DeviceType,Kokkos::Cuda>::value) {
int natom_pad = ((natom + 32 - 1) / 32) * 32; // for AoSoA layouts
bytes += natom * idxu_max * nelements * sizeof(double); // ulisttot_re
bytes += natom * idxu_max * nelements * sizeof(double); // ulisttot_im
bytes += natom_pad * idxu_max * nelements * sizeof(double) * 2; // ulisttot_pack

40
src/KSPACE/ewald.cpp
View File

@ -21,6 +21,7 @@
#include "ewald.h"
#include <mpi.h>
#include <cmath>
#include <string>
#include "atom.h"
#include "comm.h"
#include "force.h"
@ -29,6 +30,8 @@
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -88,10 +91,7 @@ Ewald::~Ewald()
void Ewald::init()
{
if (comm->me == 0) {
if (screen) fprintf(screen,"Ewald initialization ...\n");
if (logfile) fprintf(logfile,"Ewald initialization ...\n");
}
if (comm->me == 0) utils::logmesg(lmp,"Ewald initialization ...\n");
// error check
@ -185,28 +185,16 @@ void Ewald::init()
// stats
if (comm->me == 0) {
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(screen," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(logfile," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += fmt::format(" KSpace vectors: actual max1d max3d = {} {} {}\n",
kcount,kmax,kmax3d);
mesg += fmt::format(" kxmax kymax kzmax = {} {} {}\n",
kxmax,kymax,kzmax);
utils::logmesg(lmp,mesg);
}
}

40
src/KSPACE/ewald_dipole.cpp
View File

@ -18,6 +18,7 @@
#include "ewald_dipole.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include <cmath>
#include "atom.h"
#include "comm.h"
@ -29,6 +30,8 @@
#include "memory.h"
#include "error.h"
#include "update.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -63,10 +66,7 @@ EwaldDipole::~EwaldDipole()
void EwaldDipole::init()
{
if (comm->me == 0) {
if (screen) fprintf(screen,"EwaldDipole initialization ...\n");
if (logfile) fprintf(logfile,"EwaldDipole initialization ...\n");
}
if (comm->me == 0) utils::logmesg(lmp,"EwaldDipole initialization ...\n");
// error check
@ -192,28 +192,16 @@ void EwaldDipole::init()
// stats
if (comm->me == 0) {
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(screen," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(logfile," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += fmt::format(" KSpace vectors: actual max1d max3d = {} {} {}\n",
kcount,kmax,kmax3d);
mesg += fmt::format(" kxmax kymax kzmax = {} {} {}\n",
kxmax,kymax,kzmax);
utils::logmesg(lmp,mesg);
}
}

40
src/KSPACE/ewald_dipole_spin.cpp
View File

@ -18,6 +18,7 @@
#include "ewald_dipole_spin.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include <cmath>
#include "atom.h"
#include "comm.h"
@ -28,6 +29,8 @@
#include "memory.h"
#include "error.h"
#include "update.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -61,10 +64,7 @@ EwaldDipoleSpin::~EwaldDipoleSpin() {}
void EwaldDipoleSpin::init()
{
if (comm->me == 0) {
if (screen) fprintf(screen,"EwaldDipoleSpin initialization ...\n");
if (logfile) fprintf(logfile,"EwaldDipoleSpin initialization ...\n");
}
if (comm->me == 0) utils::logmesg(lmp,"EwaldDipoleSpin initialization ...\n");
// error check
@ -182,28 +182,16 @@ void EwaldDipoleSpin::init()
// stats
if (comm->me == 0) {
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(screen," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," KSpace vectors: actual max1d max3d = %d %d %d\n",
kcount,kmax,kmax3d);
fprintf(logfile," kxmax kymax kzmax = %d %d %d\n",
kxmax,kymax,kzmax);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += fmt::format(" KSpace vectors: actual max1d max3d = {} {} {}\n",
kcount,kmax,kmax3d);
mesg += fmt::format(" kxmax kymax kzmax = {} {} {}\n",
kxmax,kymax,kzmax);
utils::logmesg(lmp,mesg);
}
}

29
src/KSPACE/ewald_disp.cpp
View File

@ -18,6 +18,7 @@
#include "ewald_disp.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include <cmath>
#include "math_vector.h"
#include "math_const.h"
@ -30,6 +31,8 @@
#include "memory.h"
#include "error.h"
#include "update.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -89,10 +92,7 @@ void EwaldDisp::init()
nbox = -1;
bytes = 0.0;
if (!comm->me) {
if (screen) fprintf(screen,"EwaldDisp initialization ...\n");
if (logfile) fprintf(logfile,"EwaldDisp initialization ...\n");
}
if (!comm->me) utils::logmesg(lmp,"EwaldDisp initialization ...\n");
triclinic_check();
if (domain->dimension == 2)
@ -169,11 +169,9 @@ void EwaldDisp::init()
if (qsqsum == 0.0 && bsbsum == 0.0 && M2 == 0.0)
error->all(FLERR,"Cannot use Ewald/disp solver "
"on system with no charge, dipole, or LJ particles");
if (fabs(qsum) > SMALL && comm->me == 0) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->warning(FLERR,str);
}
if (fabs(qsum) > SMALL && comm->me == 0)
error->warning(FLERR,fmt::format("System is not charge neutral, "
"net charge = {}",qsum));
if (!function[1] && !function[2]) dispersionflag = 0;
if (!function[3]) dipoleflag = 0;
@ -229,10 +227,9 @@ void EwaldDisp::init()
}
}
if (!comm->me) {
if (screen) fprintf(screen, " G vector = %g, accuracy = %g\n", g_ewald,accuracy);
if (logfile) fprintf(logfile, " G vector = %g accuracy = %g\n", g_ewald,accuracy);
}
if (!comm->me)
utils::logmesg(lmp,fmt::format(" G vector = {}, accuracy = {}\n",
g_ewald,accuracy));
g_ewald_6 = g_ewald;
deallocate_peratom();
@ -294,10 +291,8 @@ void EwaldDisp::setup()
if (!(first_output||comm->me)) {
first_output = 1;
if (screen) fprintf(screen,
" vectors: nbox = %d, nkvec = %d\n", nbox, nkvec);
if (logfile) fprintf(logfile,
" vectors: nbox = %d, nkvec = %d\n", nbox, nkvec);
utils::logmesg(lmp,fmt::format(" vectors: nbox = {}, nkvec = {}\n",
nbox, nkvec));
}
}

31
src/KSPACE/fix_tune_kspace.cpp
View File

@ -19,6 +19,8 @@
#include <cmath>
#include <cstring>
#include <limits>
#include <string>
#include "comm.h"
#include "update.h"
#include "force.h"
#include "kspace.h"
@ -29,6 +31,8 @@
#include "neighbor.h"
#include "modify.h"
#include "compute.h"
#include "utils.h"
#include "fmt/format.h"
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
@ -146,10 +150,11 @@ void FixTuneKspace::pre_exchange()
update_kspace_style(new_kspace_style,new_acc_str);
} else if (niter == 4) {
store_old_kspace_settings();
if (screen) fprintf(screen,"ewald_time = %g\npppm_time = %g\nmsm_time = %g",
ewald_time, pppm_time, msm_time);
if (logfile) fprintf(logfile,"ewald_time = %g\npppm_time = %g\nmsm_time = %g",
ewald_time, pppm_time, msm_time);
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("ewald_time = {}\n"
"pppm_time = {}\n"
"msm_time = {}\n",
ewald_time, pppm_time, msm_time));
// switch to fastest one
strcpy(new_kspace_style,"ewald");
snprintf(new_pair_style,64,"%s/long",base_pair_style);
@ -240,8 +245,8 @@ void FixTuneKspace::update_pair_style(char *new_pair_style,
p_pair_settings_file = tmpfile();
force->pair->write_restart(p_pair_settings_file);
rewind(p_pair_settings_file);
if (screen) fprintf(screen,"Creating new pair style: %s\n",new_pair_style);
if (logfile) fprintf(logfile,"Creating new pair style: %s\n",new_pair_style);
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("Creating new pair style: {}\n",new_pair_style));
// delete old pair style and create new one
force->create_pair(new_pair_style,1);
@ -250,8 +255,9 @@ void FixTuneKspace::update_pair_style(char *new_pair_style,
double *pcutoff = (double *) force->pair->extract("cut_coul",itmp);
double current_cutoff = *pcutoff;
if (screen) fprintf(screen,"Coulomb cutoff for real space: %g\n", current_cutoff);
if (logfile) fprintf(logfile,"Coulomb cutoff for real space: %g\n", current_cutoff);
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("Coulomb cutoff for real space: {}\n",
current_cutoff));
// close temporary file
fclose(p_pair_settings_file);
@ -319,8 +325,9 @@ void FixTuneKspace::adjust_rcut(double time)
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
double current_cutoff = *p_cutoff;
if (screen) fprintf(screen,"Old Coulomb cutoff for real space: %g\n",current_cutoff);
if (logfile) fprintf(logfile,"Old Coulomb cutoff for real space: %g\n",current_cutoff);
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("Old Coulomb cutoff for real space: {}\n",
current_cutoff));
// use Brent's method from Numerical Recipes to find optimal real space cutoff
@ -390,8 +397,8 @@ void FixTuneKspace::adjust_rcut(double time)
// report the new cutoff
double *new_cutoff = (double *) force->pair->extract("cut_coul",itmp);
current_cutoff = *new_cutoff;
if (screen) fprintf(screen,"Adjusted Coulomb cutoff for real space: %g\n", current_cutoff);
if (logfile) fprintf(logfile,"Adjusted Coulomb cutoff for real space: %g\n", current_cutoff);
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("Adjusted Coulomb cutoff for real space: {}\n", current_cutoff));
store_old_kspace_settings();
update_pair_style(new_pair_style,pair_cut_coul);

22
src/KSPACE/msm_cg.cpp
View File

@ -26,6 +26,8 @@
#include "force.h"
#include "neighbor.h"
#include "memory.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -140,20 +142,12 @@ void MSMCG::compute(int eflag, int vflag)
charged_frac = static_cast<double>(charged_all) * 100.0
/ static_cast<double>(atom->natoms);
if (me == 0) {
if (screen)
fprintf(screen,
" MSM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
if (logfile)
fprintf(logfile,
" MSM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
}
if (me == 0)
utils::logmesg(lmp,fmt::format(" MSM/cg optimization cutoff: {}\n"
" Total charged atoms: {:.1f}%\n"
" Min/max charged atoms/proc: {:.1f}%"
" {:.1f}%\n",smallq,
charged_frac,charged_fmin,charged_fmax));
}
// only need to rebuild this list after a neighbor list update

68
src/KSPACE/pppm.cpp
View File

@ -21,6 +21,7 @@
#include "pppm.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include <cmath>
#include "atom.h"
#include "comm.h"
@ -35,6 +36,8 @@
#include "remap_wrap.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
#include "math_const.h"
#include "math_special.h"
@ -185,10 +188,7 @@ PPPM::~PPPM()
void PPPM::init()
{
if (me == 0) {
if (screen) fprintf(screen,"PPPM initialization ...\n");
if (logfile) fprintf(logfile,"PPPM initialization ...\n");
}
if (me == 0) utils::logmesg(lmp,"PPPM initialization ...\n");
// error check
@ -203,13 +203,14 @@ void PPPM::init()
if (domain->triclinic && slabflag)
error->all(FLERR,"Cannot (yet) use PPPM with triclinic box and "
"slab correction");
if (domain->dimension == 2) error->all(FLERR,
"Cannot use PPPM with 2d simulation");
if (domain->dimension == 2)
error->all(FLERR,"Cannot use PPPM with 2d simulation");
if (comm->style != 0)
error->universe_all(FLERR,"PPPM can only currently be used with "
"comm_style brick");
if (!atom->q_flag) error->all(FLERR,"Kspace style requires atom attribute q");
if (!atom->q_flag)
error->all(FLERR,"Kspace style requires atom attribute q");
if (slabflag == 0 && domain->nonperiodic > 0)
error->all(FLERR,"Cannot use non-periodic boundaries with PPPM");
@ -219,11 +220,8 @@ void PPPM::init()
error->all(FLERR,"Incorrect boundaries with slab PPPM");
}
if (order < 2 || order > MAXORDER) {
char str[128];
sprintf(str,"PPPM order cannot be < 2 or > than %d",MAXORDER);
error->all(FLERR,str);
}
if (order < 2 || order > MAXORDER)
error->all(FLERR,fmt::format("PPPM order cannot be < 2 or > {}",MAXORDER));
// compute two charge force
@ -246,10 +244,7 @@ void PPPM::init()
qdist = 0.0;
if (tip4pflag) {
if (me == 0) {
if (screen) fprintf(screen," extracting TIP4P info from pair style\n");
if (logfile) fprintf(logfile," extracting TIP4P info from pair style\n");
}
if (me == 0) utils::logmesg(lmp," extracting TIP4P info from pair style\n");
double *p_qdist = (double *) force->pair->extract("qdist",itmp);
int *p_typeO = (int *) force->pair->extract("typeO",itmp);
@ -351,31 +346,17 @@ void PPPM::init()
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(screen," stencil order = %d\n",order);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n");
fprintf(screen," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(logfile," stencil order = %d\n",order);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n");
fprintf(logfile," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" grid = {} {} {}\n",nx_pppm,ny_pppm,nz_pppm);
mesg += fmt::format(" stencil order = {}\n",order);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += " using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n";
mesg += fmt::format(" 3d grid and FFT values/proc = {} {}\n",
ngrid_max,nfft_both_max);
utils::logmesg(lmp,mesg);
}
// allocate K-space dependent memory
@ -1291,10 +1272,7 @@ void PPPM::adjust_gewald()
g_ewald -= dx;
if (fabs(newton_raphson_f()) < SMALL) return;
}
char str[128];
sprintf(str, "Could not compute g_ewald");
error->all(FLERR, str);
error->all(FLERR, "Could not compute g_ewald");
}
/* ----------------------------------------------------------------------

22
src/KSPACE/pppm_cg.cpp
View File

@ -28,6 +28,8 @@
#include "memory.h"
#include "math_const.h"
#include "remap.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -151,20 +153,12 @@ void PPPMCG::compute(int eflag, int vflag)
charged_frac = static_cast<double>(charged_all) * 100.0
/ static_cast<double>(atom->natoms);
if (me == 0) {
if (screen)
fprintf(screen,
" PPPM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
if (logfile)
fprintf(logfile,
" PPPM/cg optimization cutoff: %g\n"
" Total charged atoms: %.1f%%\n"
" Min/max charged atoms/proc: %.1f%% %.1f%%\n",
smallq,charged_frac,charged_fmin,charged_fmax);
}
if (me == 0)
utils::logmesg(lmp,fmt::format(" PPPM/cg optimization cutoff: {}\n"
" Total charged atoms: {:.1f}%\n"
" Min/max charged atoms/proc: {:.1f}%"
" {:.1f}%\n",smallq,
charged_frac,charged_fmin,charged_fmax));
}
// only need to rebuild this list after a neighbor list update

58
src/KSPACE/pppm_dipole.cpp
View File

@ -18,6 +18,7 @@
#include "pppm_dipole.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include <cmath>
#include "atom.h"
#include "comm.h"
@ -30,6 +31,8 @@
#include "memory.h"
#include "error.h"
#include "update.h"
#include "utils.h"
#include "fmt/format.h"
#include "math_const.h"
#include "math_special.h"
@ -102,10 +105,7 @@ PPPMDipole::~PPPMDipole()
void PPPMDipole::init()
{
if (me == 0) {
if (screen) fprintf(screen,"PPPMDipole initialization ...\n");
if (logfile) fprintf(logfile,"PPPMDipole initialization ...\n");
}
if (me == 0) utils::logmesg(lmp,"PPPMDipole initialization ...\n");
// error check
@ -143,11 +143,9 @@ void PPPMDipole::init()
error->all(FLERR,"Incorrect boundaries with slab PPPMDipole");
}
if (order < 2 || order > MAXORDER) {
char str[128];
sprintf(str,"PPPMDipole order cannot be < 2 or > than %d",MAXORDER);
error->all(FLERR,str);
}
if (order < 2 || order > MAXORDER)
error->all(FLERR,fmt::format("PPPMDipole order cannot be < 2 or > {}",
MAXORDER));
// compute two charge force
@ -246,37 +244,17 @@ void PPPMDipole::init()
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
#ifdef FFT_SINGLE
const char fft_prec[] = "single";
#else
const char fft_prec[] = "double";
#endif
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(screen," stencil order = %d\n",order);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," using %s precision FFTs\n",fft_prec);
fprintf(screen," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(logfile," stencil order = %d\n",order);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," using %s precision FFTs\n",fft_prec);
fprintf(logfile," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" grid = {} {} {}\n",nx_pppm,ny_pppm,nz_pppm);
mesg += fmt::format(" stencil order = {}\n",order);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += " using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n";
mesg += fmt::format(" 3d grid and FFT values/proc = {} {}\n",
ngrid_max,nfft_both_max);
utils::logmesg(lmp,mesg);
}
// allocate K-space dependent memory

58
src/KSPACE/pppm_dipole_spin.cpp
View File

@ -18,6 +18,7 @@
#include "pppm_dipole_spin.h"
#include <mpi.h>
#include <cstring>
#include <string>
#include "atom.h"
#include "comm.h"
#include "gridcomm.h"
@ -27,6 +28,8 @@
#include "memory.h"
#include "error.h"
#include "update.h"
#include "utils.h"
#include "fmt/format.h"
#include "math_const.h"
@ -87,10 +90,7 @@ PPPMDipoleSpin::~PPPMDipoleSpin()
void PPPMDipoleSpin::init()
{
if (me == 0) {
if (screen) fprintf(screen,"PPPMDipoleSpin initialization ...\n");
if (logfile) fprintf(logfile,"PPPMDipoleSpin initialization ...\n");
}
if (me == 0) utils::logmesg(lmp,"PPPMDipoleSpin initialization ...\n");
// error check
@ -123,11 +123,9 @@ void PPPMDipoleSpin::init()
error->all(FLERR,"Incorrect boundaries with slab PPPMDipoleSpin");
}
if (order < 2 || order > MAXORDER) {
char str[128];
sprintf(str,"PPPMDipoleSpin order cannot be < 2 or > than %d",MAXORDER);
error->all(FLERR,str);
}
if (order < 2 || order > MAXORDER)
error->all(FLERR,fmt::format("PPPMDipoleSpin order cannot be < 2 or > {}",
MAXORDER));
// compute two charge force
@ -226,37 +224,17 @@ void PPPMDipoleSpin::init()
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
#ifdef FFT_SINGLE
const char fft_prec[] = "single";
#else
const char fft_prec[] = "double";
#endif
if (screen) {
fprintf(screen," G vector (1/distance) = %g\n",g_ewald);
fprintf(screen," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(screen," stencil order = %d\n",order);
fprintf(screen," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(screen," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(screen," using %s precision FFTs\n",fft_prec);
fprintf(screen," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
if (logfile) {
fprintf(logfile," G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(logfile," stencil order = %d\n",order);
fprintf(logfile," estimated absolute RMS force accuracy = %g\n",
estimated_accuracy);
fprintf(logfile," estimated relative force accuracy = %g\n",
estimated_accuracy/two_charge_force);
fprintf(logfile," using %s precision FFTs\n",fft_prec);
fprintf(logfile," 3d grid and FFT values/proc = %d %d\n",
ngrid_max,nfft_both_max);
}
std::string mesg = fmt::format(" G vector (1/distance) = {}\n",g_ewald);
mesg += fmt::format(" grid = {} {} {}\n",nx_pppm,ny_pppm,nz_pppm);
mesg += fmt::format(" stencil order = {}\n",order);
mesg += fmt::format(" estimated absolute RMS force accuracy = {}\n",
estimated_accuracy);
mesg += fmt::format(" estimated relative force accuracy = {}\n",
estimated_accuracy/two_charge_force);
mesg += " using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n";
mesg += fmt::format(" 3d grid and FFT values/proc = {} {}\n",
ngrid_max,nfft_both_max);
utils::logmesg(lmp,mesg);
}
// allocate K-space dependent memory

168
src/KSPACE/pppm_disp.cpp
View File

@ -20,6 +20,7 @@
#include <mpi.h>
#include <cstring>
#include <cmath>
#include <string>
#include "math_const.h"
#include "atom.h"
#include "comm.h"
@ -34,6 +35,8 @@
#include "remap_wrap.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -254,10 +257,7 @@ PPPMDisp::~PPPMDisp()
void PPPMDisp::init()
{
if (me == 0) {
if (screen) fprintf(screen,"PPPMDisp initialization ...\n");
if (logfile) fprintf(logfile,"PPPMDisp initialization ...\n");
}
if (me == 0) utils::logmesg(lmp,"PPPMDisp initialization ...\n");
// error check
@ -277,11 +277,9 @@ void PPPMDisp::init()
error->all(FLERR,"Incorrect boundaries with slab PPPMDisp");
}
if (order > MAXORDER || order_6 > MAXORDER) {
char str[128];
sprintf(str,"PPPMDisp coulomb order cannot be greater than %d",MAXORDER);
error->all(FLERR,str);
}
if (order > MAXORDER || order_6 > MAXORDER)
error->all(FLERR,fmt::format("PPPMDisp coulomb or dispersion order cannot"
" be greater than {}",MAXORDER));
// compute two charge force
@ -318,7 +316,6 @@ void PPPMDisp::init()
for (int i=0; i<=EWALD_MAXORDER; ++i) // transcribe order
if (ewald_order&(1<<i)) { // from pair_style
int k=0;
char str[128];
switch (i) {
case 1:
k = 0; break;
@ -330,9 +327,9 @@ void PPPMDisp::init()
else error->all(FLERR,"Unsupported mixing rule in kspace_style pppm/disp");
break;
default:
sprintf(str, "Unsupported order in kspace_style "
"pppm/disp, pair_style %s", force->pair_style);
error->all(FLERR,str);
error->all(FLERR,std::string("Unsupported order in kspace_style "
"pppm/disp, pair_style ")
+ force->pair_style);
}
function[k] = 1;
}
@ -340,11 +337,8 @@ void PPPMDisp::init()
// warn, if function[0] is not set but charge attribute is set!
if (!function[0] && atom->q_flag && me == 0) {
char str[128];
sprintf(str, "Charges are set, but coulombic solver is not used");
error->warning(FLERR, str);
}
if (!function[0] && atom->q_flag && me == 0)
error->warning(FLERR, "Charges are set, but coulombic solver is not used");
// show error message if pppm/disp is not used correctly
@ -481,31 +475,19 @@ void PPPMDisp::init()
MPI_Allreduce(&nfft_both,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (screen) {
fprintf(screen," Coulomb G vector (1/distance)= %g\n",g_ewald);
fprintf(screen," Coulomb grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(screen," Coulomb stencil order = %d\n",order);
fprintf(screen," Coulomb estimated absolute RMS force accuracy = %g\n",
acc);
fprintf(screen," Coulomb estimated relative force accuracy = %g\n",
acc/two_charge_force);
fprintf(screen," using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n");
fprintf(screen," 3d grid and FFT values/proc = %d %d\n",
ngrid_max, nfft_both_max);
}
if (logfile) {
fprintf(logfile," Coulomb G vector (1/distance) = %g\n",g_ewald);
fprintf(logfile," Coulomb grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
fprintf(logfile," Coulomb stencil order = %d\n",order);
fprintf(logfile,
" Coulomb estimated absolute RMS force accuracy = %g\n",
acc);
fprintf(logfile," Coulomb estimated relative force accuracy = %g\n",
acc/two_charge_force);
fprintf(logfile," using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n");
fprintf(logfile," 3d grid and FFT values/proc = %d %d\n",
ngrid_max, nfft_both_max);
}
std::string mesg = fmt::format(" Coulomb G vector (1/distance)= {}\n",
g_ewald);
mesg += fmt::format(" Coulomb grid = {} {} {}\n",
nx_pppm,ny_pppm,nz_pppm);
mesg += fmt::format(" Coulomb stencil order = {}\n",order);
mesg += fmt::format(" Coulomb estimated absolute RMS force accuracy "
"= {}\n",acc);
mesg += fmt::format(" Coulomb estimated relative force accuracy = {}\n",
acc/two_charge_force);
mesg += " using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n";
mesg += fmt::format(" 3d grid and FFT values/proc = {} {}\n",
ngrid_max, nfft_both_max);
utils::logmesg(lmp,mesg);
}
}
@ -573,46 +555,19 @@ void PPPMDisp::init()
MPI_Allreduce(&nfft_both_6,&nfft_both_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
#ifdef FFT_SINGLE
const char fft_prec[] = "single";
#else
const char fft_prec[] = "double";
#endif
if (screen) {
fprintf(screen," Dispersion G vector (1/distance)= %g\n",g_ewald_6);
fprintf(screen," Dispersion grid = %d %d %d\n",
nx_pppm_6,ny_pppm_6,nz_pppm_6);
fprintf(screen," Dispersion stencil order = %d\n",order_6);
fprintf(screen," Dispersion estimated absolute "
"RMS force accuracy = %g\n",acc);
fprintf(screen," Dispersion estimated absolute "
"real space RMS force accuracy = %g\n",acc_real);
fprintf(screen," Dispersion estimated absolute "
"kspace RMS force accuracy = %g\n",acc_kspace);
fprintf(screen," Dispersion estimated relative force accuracy = %g\n",
acc/two_charge_force);
fprintf(screen," using %s precision FFTs\n",fft_prec);
fprintf(screen," 3d grid and FFT values/proc dispersion = %d %d\n",
ngrid_max,nfft_both_max);
}
if (logfile) {
fprintf(logfile," Dispersion G vector (1/distance) = %g\n",g_ewald_6);
fprintf(logfile," Dispersion grid = %d %d %d\n",
nx_pppm_6,ny_pppm_6,nz_pppm_6);
fprintf(logfile," Dispersion stencil order = %d\n",order_6);
fprintf(logfile," Dispersion estimated absolute "
"RMS force accuracy = %g\n",acc);
fprintf(logfile," Dispersion estimated absolute "
"real space RMS force accuracy = %g\n",acc_real);
fprintf(logfile," Dispersion estimated absolute "
"kspace RMS force accuracy = %g\n",acc_kspace);
fprintf(logfile," Disperion estimated relative force accuracy = %g\n",
acc/two_charge_force);
fprintf(logfile," using %s precision FFTs\n",fft_prec);
fprintf(logfile," 3d grid and FFT values/proc dispersion = %d %d\n",
ngrid_max,nfft_both_max);
}
std::string mesg = fmt::format(" Dispersion G vector (1/distance)= "
"{}\n", g_ewald_6);
mesg += fmt::format(" Dispersion grid = {} {} {}\n",
nx_pppm_6,ny_pppm_6,nz_pppm_6);
mesg += fmt::format(" Dispersion stencil order = {}\n",order_6);
mesg += fmt::format(" Dispersion estimated absolute RMS force accuracy "
"= {}\n",acc);
mesg += fmt::format(" Dispersion estimated relative force accuracy "
"= {}\n",acc/two_charge_force);
mesg += " using " LMP_FFT_PREC " precision " LMP_FFT_LIB "\n";
mesg += fmt::format(" 3d grid and FFT values/proc = {} {}\n",
ngrid_max, nfft_both_max);
utils::logmesg(lmp,mesg);
}
}
@ -1285,10 +1240,8 @@ void PPPMDisp::init_coeffs() // local pair coeffs
delete [] B;
B = NULL;
if (function[3] + function[2]) { // no mixing rule or arithmetic
if (function[2] && me == 0) {
if (screen) fprintf(screen," Optimizing splitting of Dispersion coefficients\n");
if (logfile) fprintf(logfile," Optimizing splitting of Dispersion coefficients\n");
}
if (function[2] && me == 0)
utils::logmesg(lmp," Optimizing splitting of Dispersion coefficients\n");
// allocate data for eigenvalue decomposition
double **A=NULL;
@ -1349,11 +1302,9 @@ void PPPMDisp::init_coeffs() // local pair coeffs
}
err = bmax/amax;
if (err > 1.0e-4) {
char str[128];
sprintf(str,"Estimated error in splitting of dispersion coeffs is %g",err);
error->warning(FLERR, str);
}
if (err > 1.0e-4 && comm->me == 0)
error->warning(FLERR,fmt::format("Estimated error in splitting of "
"dispersion coeffs is {}",err));
// set B
B = new double[nsplit*n+nsplit];
for (int i = 0; i< nsplit; i++) {
@ -1374,31 +1325,24 @@ void PPPMDisp::init_coeffs() // local pair coeffs
function[3] = 0;
function[2] = 0;
function[1] = 1;
if (me == 0) {
if (screen) fprintf(screen," Using geometric mixing for reciprocal space\n");
if (logfile) fprintf(logfile," Using geometric mixing for reciprocal space\n");
}
if (me == 0)
utils::logmesg(lmp," Using geometric mixing for reciprocal space\n");
}
if (function[2] && nsplit <= 6) {
if (me == 0) {
if (screen) fprintf(screen," Using %d instead of 7 structure factors\n",nsplit);
if (logfile) fprintf(logfile," Using %d instead of 7 structure factors\n",nsplit);
}
if (me == 0)
utils::logmesg(lmp,fmt::format(" Using {} instead of 7 structure "
"factors\n",nsplit));
function[3] = 1;
function[2] = 0;
}
if (function[2] && (nsplit > 6)) {
if (me == 0) {
if (screen) fprintf(screen," Using 7 structure factors\n");
if (logfile) fprintf(logfile," Using 7 structure factors\n");
}
if (me == 0) utils::logmesg(lmp," Using 7 structure factors\n");
if ( B ) delete [] B;
}
if (function[3]) {
if (me == 0) {
if (screen) fprintf(screen," Using %d structure factors\n",nsplit);
if (logfile) fprintf(logfile," Using %d structure factors\n",nsplit);
}
if (me == 0)
utils::logmesg(lmp,fmt::format(" Using {} structure factors\n",
nsplit));
if (nsplit > 9) error->warning(FLERR, "Simulations might be very slow because of large number of structure factors");
}
@ -2884,9 +2828,7 @@ void PPPMDisp::adjust_gewald()
// Failed to converge
char str[128];
sprintf(str, "Could not compute g_ewald");
error->all(FLERR, str);
error->all(FLERR, "Could not compute g_ewald");
}
@ -3517,9 +3459,7 @@ void PPPMDisp::adjust_gewald_6()
// Failed to converge
char str[128];
sprintf(str, "Could not adjust g_ewald_6");
error->all(FLERR, str);
error->all(FLERR, "Could not adjust g_ewald_6");
}

24
src/MANYBODY/fix_qeq_comb.cpp
View File

@ -31,6 +31,7 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -64,11 +65,9 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
if (iarg+2 > narg) error->all(FLERR,"Illegal fix qeq/comb command");
if (me == 0) {
fp = fopen(arg[iarg+1],"w");
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
error->one(FLERR,str);
}
if (fp == NULL)
error->one(FLERR,std::string("Cannot open fix qeq/comb file ")
+ arg[iarg+1]);
}
iarg += 2;
} else error->all(FLERR,"Illegal fix qeq/comb command");
@ -191,8 +190,7 @@ void FixQEQComb::post_force(int /*vflag*/)
// charge-equilibration loop
if (me == 0 && fp)
fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
update->ntimestep);
fmt::print(fp,"Charge equilibration on step {}\n", update->ntimestep);
heatpq = 0.05;
qmass = 0.016;
@ -253,9 +251,9 @@ void FixQEQComb::post_force(int /*vflag*/)
if (enegchk <= precision && enegmax <= 100.0*precision) break;
if (me == 0 && fp)
fprintf(fp," iteration: %d, enegtot %.6g, "
"enegmax %.6g, fq deviation: %.6g\n",
iloop,enegtot,enegmax,enegchk);
fmt::print(fp," iteration: {}, enegtot {:.6g}, "
"enegmax {:.6g}, fq deviation: {:.6g}\n",
iloop,enegtot,enegmax,enegchk);
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
@ -266,10 +264,10 @@ void FixQEQComb::post_force(int /*vflag*/)
if (me == 0 && fp) {
if (iloop == loopmax)
fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
fmt::print(fp,"Charges did not converge in {} iterations\n",iloop);
else
fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
iloop,enegchk);
fmt::print(fp,"Charges converged in {} iterations to {:.10f} tolerance\n",
iloop,enegchk);
}
}

10
src/MANYBODY/pair_bop.cpp
View File

@ -48,6 +48,7 @@
#include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -689,12 +690,9 @@ void PairBOP::init_style()
// check that user sets comm->cutghostuser to 3x the max BOP cutoff
if (comm->cutghostuser < 3.0*cutmax - EPSILON) {
char str[128];
sprintf(str,"Pair style bop requires comm ghost cutoff "
"at least 3x larger than %g",cutmax);
error->all(FLERR,str);
}
if (comm->cutghostuser < 3.0*cutmax - EPSILON)
error->all(FLERR,fmt::format("Pair style bop requires comm ghost cutoff "
"at least 3x larger than {}",cutmax));
// need a full neighbor list and neighbors of ghosts

14
src/MANYBODY/pair_comb.cpp
View File

@ -493,15 +493,17 @@ void PairComb::coeff(int narg, char **arg)
// generate streitz-mintmire direct 1/r energy look-up table
if (comm->me == 0 && screen) fprintf(screen,"Pair COMB:\n");
if (comm->me == 0 && screen)
fprintf(screen," generating Coulomb integral lookup table ...\n");
fputs("Pair COMB:\n"
" generating Coulomb integral lookup table ...\n",screen);
sm_table();
if (cor_flag && comm->me == 0 && screen)
fprintf(screen," will apply over-coordination correction ...\n");
if (!cor_flag&& comm->me == 0 && screen)
fprintf(screen," will not apply over-coordination correction ...\n");
if (comm->me == 0 && screen) {
if (cor_flag)
fputs(" will apply over-coordination correction ...\n",screen);
else
fputs(" will not apply over-coordination correction ...\n",screen);
}
// clear setflag since coeff() called once with I,J = * *

23
src/MANYBODY/pair_comb3.cpp
View File

@ -36,6 +36,7 @@
#include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace MathConst;
@ -165,17 +166,13 @@ void PairComb3::settings(int narg, char **arg)
if (narg != 1) error->all(FLERR,"Illegal pair_style command");
if (strcmp(arg[0],"polar_on") == 0) {
pol_flag = 1;
if (comm->me == 0 && screen) fprintf(screen,
" PairComb3: polarization is on \n");
} else if (strcmp(arg[0],"polar_off") == 0) {
pol_flag = 0;
if (comm->me == 0 && screen) fprintf(screen,
" PairComb3: polarization is off \n");
} else {
error->all(FLERR,"Illegal pair_style command");
}
if (strcmp(arg[0],"polar_on") == 0) pol_flag = 1;
else if (strcmp(arg[0],"polar_off") == 0) pol_flag = 0;
else error->all(FLERR,"Illegal pair_style command");
if (comm->me == 0 && screen)
fmt::print(screen," PairComb3: polarization is {} \n",
pol_flag ? "on" : "off");
}
/* ----------------------------------------------------------------------
@ -210,8 +207,8 @@ void PairComb3::coeff(int narg, char **arg)
nelements = 0;
for (i = 3; i < narg; i++) {
if ((strcmp(arg[i],"C") == 0) && (cflag == 0)) {
if (comm->me == 0 && screen) fprintf(screen,
" PairComb3: Found C: reading additional library file\n");
if (comm->me == 0 && screen)
fputs(" PairComb3: Found C: reading additional library file\n",screen);
read_lib();
cflag = 1;
}

17
src/MANYBODY/pair_eam.cpp
View File

@ -43,6 +43,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
manybody_flag = 1;
embedstep = -1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nmax = 0;
rho = NULL;
@ -466,8 +467,13 @@ void PairEAM::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(lmp, filename, "EAM", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
@ -492,6 +498,14 @@ void PairEAM::read_file(char *filename)
reader.next_dvector(&file->frho[1], file->nrho);
reader.next_dvector(&file->zr[1], file->nr);
reader.next_dvector(&file->rhor[1], file->nr);
if (unit_convert) {
const double sqrt_conv = sqrt(conversion_factor);
for (int i = 1; i <= file->nrho; ++i)
file->frho[i] *= conversion_factor;
for (int j = 1; j <= file->nr; ++j)
file->zr[j] *= sqrt_conv;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}
@ -655,6 +669,7 @@ void PairEAM::file2array()
// create a z2r array for each file against other files, only for I >= J
// interpolate zri and zrj to a single grid and cutoff
// final z2r includes unit conversion of 27.2 eV/Hartree and 0.529 Ang/Bohr
double zri,zrj;

15
src/MANYBODY/pair_eam_alloy.cpp
View File

@ -119,8 +119,13 @@ void PairEAMAlloy::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(lmp, filename, "EAMAlloy", unit_convert_flag);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -165,11 +170,19 @@ void PairEAMAlloy::read_file(char *filename)
reader.next_dvector(&file->frho[i][1], file->nrho);
reader.next_dvector(&file->rhor[i][1], file->nr);
if (unit_convert) {
for (int j = 1; j < file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
}
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k < file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

14
src/MANYBODY/pair_eam_cd.cpp
View File

@ -29,6 +29,8 @@
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -40,6 +42,7 @@ PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion)
{
single_enable = 0;
restartinfo = 0;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
rhoB = NULL;
D_values = NULL;
@ -500,12 +503,11 @@ void PairEAMCD::read_h_coeff(char *filename)
FILE *fptr;
char line[MAXLINE];
char nextline[MAXLINE];
fptr = force->open_potential(filename);
if (fptr == NULL) {
char str[128];
snprintf(str,128,"Cannot open EAM potential file %s", filename);
error->one(FLERR,str);
}
int convert_flag = unit_convert_flag;
fptr = force->open_potential(filename, &convert_flag);
if (fptr == NULL)
error->one(FLERR,fmt::format("Cannot open EAMCD potential file {}",
filename));
// h coefficients are stored at the end of the file.
// Skip to last line of file.

17
src/MANYBODY/pair_eam_fs.cpp
View File

@ -119,8 +119,15 @@ void PairEAMFS::read_file(char *filename)
// read potential file
if(comm->me == 0) {
PotentialFileReader reader(lmp, filename, "EAM");
PotentialFileReader reader(lmp, filename, "EAMFS", unit_convert_flag);
// transparently convert units for supported conversions
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
try {
reader.skip_line();
reader.skip_line();
@ -164,6 +171,10 @@ void PairEAMFS::read_file(char *filename)
file->mass[i] = values.next_double();
reader.next_dvector(&file->frho[i][1], file->nrho);
if (unit_convert) {
for (int j = 1; j <= file->nrho; ++j)
file->frho[i][j] *= conversion_factor;
}
for (int j = 0; j < file->nelements; j++) {
reader.next_dvector(&file->rhor[i][j][1], file->nr);
@ -173,6 +184,10 @@ void PairEAMFS::read_file(char *filename)
for (int i = 0; i < file->nelements; i++) {
for (int j = 0; j <= i; j++) {
reader.next_dvector(&file->z2r[i][j][1], file->nr);
if (unit_convert) {
for (int k = 1; k <= file->nr; ++k)
file->z2r[i][j][k] *= conversion_factor;
}
}
}
} catch (TokenizerException & e) {

43
src/MANYBODY/pair_eim.cpp
View File

@ -41,6 +41,7 @@ PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
setfl = NULL;
nmax = 0;
@ -477,7 +478,7 @@ void PairEIM::read_file(char *filename)
// read potential file
if( comm->me == 0) {
EIMPotentialFileReader reader(lmp, filename);
EIMPotentialFileReader reader(lmp, filename, unit_convert_flag);
reader.get_global(setfl);
@ -1050,14 +1051,18 @@ double PairEIM::memory_usage()
return bytes;
}
EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS * lmp, const std::string & filename) :
EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS *lmp,
const std::string &filename,
const int auto_convert) :
Pointers(lmp), filename(filename)
{
if (comm->me != 0) {
error->one(FLERR, "EIMPotentialFileReader should only be called by proc 0!");
}
FILE * fp = force->open_potential(filename.c_str());
int unit_convert = auto_convert;
FILE *fp = force->open_potential(filename.c_str(), &unit_convert);
conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert);
if (fp == NULL) {
error->one(FLERR, fmt::format("cannot open EIM potential file {}", filename));
@ -1186,7 +1191,7 @@ void EIMPotentialFileReader::parse(FILE * fp)
PairData data;
data.rcutphiA = values.next_double();
data.rcutphiR = values.next_double();
data.Eb = values.next_double();
data.Eb = values.next_double() * conversion_factor;
data.r0 = values.next_double();
data.alpha = values.next_double();
data.beta = values.next_double();
@ -1194,7 +1199,7 @@ void EIMPotentialFileReader::parse(FILE * fp)
data.Asigma = values.next_double();
data.rq = values.next_double();
data.rcutsigma = values.next_double();
data.Ac = values.next_double();
data.Ac = values.next_double() * conversion_factor;
data.zeta = values.next_double();
data.rs = values.next_double();
@ -1217,20 +1222,18 @@ void EIMPotentialFileReader::parse(FILE * fp)
}
}
void EIMPotentialFileReader::get_global(PairEIM::Setfl * setfl) {
void EIMPotentialFileReader::get_global(PairEIM::Setfl *setfl) {
setfl->division = division;
setfl->rbig = rbig;
setfl->rsmall = rsmall;
}
void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const std::string & name) {
if (elements.find(name) == elements.end()) {
char str[128];
snprintf(str, 128, "Element %s not defined in EIM potential file", name.c_str());
error->one(FLERR, str);
}
void EIMPotentialFileReader::get_element(PairEIM::Setfl *setfl, int i,
const std::string &name) {
if (elements.find(name) == elements.end())
error->one(FLERR,"Element " + name + " not defined in EIM potential file");
ElementData & data = elements[name];
ElementData &data = elements[name];
setfl->ielement[i] = data.ielement;
setfl->mass[i] = data.mass;
setfl->negativity[i] = data.negativity;
@ -1240,16 +1243,16 @@ void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const st
setfl->q0[i] = data.q0;
}
void EIMPotentialFileReader::get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB) {
void EIMPotentialFileReader::get_pair(PairEIM::Setfl *setfl, int ij,
const std::string &elemA,
const std::string &elemB) {
auto p = get_pair(elemA, elemB);
if (pairs.find(p) == pairs.end()) {
char str[128];
snprintf(str, 128, "Pair (%s, %s) not defined in EIM potential file", elemA.c_str(), elemB.c_str());
error->one(FLERR, str);
}
if (pairs.find(p) == pairs.end())
error->one(FLERR,"Element pair (" + elemA + ", " + elemB
+ ") is not defined in EIM potential file");
PairData & data = pairs[p];
PairData &data = pairs[p];
setfl->rcutphiA[ij] = data.rcutphiA;
setfl->rcutphiR[ij] = data.rcutphiR;
setfl->Eb[ij] = data.Eb;

18
src/MANYBODY/pair_eim.h
View File

@ -98,17 +98,21 @@ class EIMPotentialFileReader : protected Pointers {
std::string filename;
static const int MAXLINE = 1024;
char line[MAXLINE];
double conversion_factor;
void parse(FILE * fp);
char * next_line(FILE * fp);
std::pair<std::string, std::string> get_pair(const std::string & a, const std::string & b);
void parse(FILE *fp);
char *next_line(FILE *fp);
std::pair<std::string, std::string> get_pair(const std::string &a,
const std::string &b);
public:
EIMPotentialFileReader(class LAMMPS* lmp, const std::string & filename);
EIMPotentialFileReader(class LAMMPS* lmp, const std::string &filename,
const int auto_convert=0);
void get_global(PairEIM::Setfl * setfl);
void get_element(PairEIM::Setfl * setfl, int i, const std::string & name);
void get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB);
void get_global(PairEIM::Setfl *setfl);
void get_element(PairEIM::Setfl *setfl, int i, const std::string &name);
void get_pair(PairEIM::Setfl *setfl, int ij,
const std::string &elemA, const std::string &elemB);
private:
// potential parameters

13
src/MANYBODY/pair_gw.cpp
View File

@ -48,6 +48,7 @@ PairGW::PairGW(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nelements = 0;
elements = NULL;
@ -375,9 +376,14 @@ void PairGW::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "GW");
PotentialFileReader reader(lmp, file, "GW", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -427,6 +433,11 @@ void PairGW::read_file(char *file)
params[nparams].lam1 = values.next_double();
params[nparams].biga = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

12
src/MANYBODY/pair_gw_zbl.cpp
View File

@ -70,9 +70,14 @@ void PairGWZBL::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "GW/ZBL");
PotentialFileReader reader(lmp, file, "GW/ZBL", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -126,6 +131,11 @@ void PairGWZBL::read_file(char *file)
params[nparams].ZBLcut = values.next_double();
params[nparams].ZBLexpscale = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

10
src/MANYBODY/pair_nb3b_harmonic.cpp
View File

@ -47,6 +47,7 @@ PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nelements = 0;
elements = NULL;
@ -291,9 +292,14 @@ void PairNb3bHarmonic::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "nb3b/harmonic");
PotentialFileReader reader(lmp, file, "nb3b/harmonic", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -331,6 +337,8 @@ void PairNb3bHarmonic::read_file(char *file)
params[nparams].k_theta = values.next_double();
params[nparams].theta0 = values.next_double();
params[nparams].cutoff = values.next_double();
if (unit_convert) params[nparams].k_theta *= conversion_factor;
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

14
src/MANYBODY/pair_sw.cpp
View File

@ -44,6 +44,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nelements = 0;
elements = NULL;
@ -355,9 +356,16 @@ void PairSW::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Stillinger-Weber");
PotentialFileReader reader(lmp, file, "Stillinger-Weber",
unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -407,6 +415,10 @@ void PairSW::read_file(char *file)
error->one(FLERR, e.what());
}
if (unit_convert) {
params[nparams].epsilon *= conversion_factor;
}
if (params[nparams].epsilon < 0.0 || params[nparams].sigma < 0.0 ||
params[nparams].littlea < 0.0 || params[nparams].lambda < 0.0 ||
params[nparams].gamma < 0.0 || params[nparams].biga < 0.0 ||

16
src/MANYBODY/pair_tersoff.cpp
View File

@ -49,6 +49,7 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nelements = 0;
elements = NULL;
@ -400,9 +401,14 @@ void PairTersoff::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Tersoff");
char * line;
PotentialFileReader reader(lmp, file, "Tersoff", unit_convert_flag);
char *line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -426,7 +432,6 @@ void PairTersoff::read_file(char *file)
if (kname == elements[kelement]) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
@ -453,6 +458,11 @@ void PairTersoff::read_file(char *file)
params[nparams].lam1 = values.next_double();
params[nparams].biga = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

12
src/MANYBODY/pair_tersoff_mod.cpp
View File

@ -53,9 +53,14 @@ void PairTersoffMOD::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Tersoff");
PotentialFileReader reader(lmp, file, "TersoffMod", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -109,6 +114,11 @@ void PairTersoffMOD::read_file(char *file)
params[nparams].c4 = values.next_double();
params[nparams].c5 = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

13
src/MANYBODY/pair_tersoff_mod_c.cpp
View File

@ -44,9 +44,14 @@ void PairTersoffMODC::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Tersoff");
PotentialFileReader reader(lmp, file, "TersoffModC", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -101,6 +106,12 @@ void PairTersoffMODC::read_file(char *file)
params[nparams].c5 = values.next_double();
params[nparams].c0 = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
params[nparams].c0 *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

14
src/MANYBODY/pair_tersoff_zbl.cpp
View File

@ -71,9 +71,15 @@ void PairTersoffZBL::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Tersoff");
PotentialFileReader reader(lmp, file, "TersoffZBL", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -97,7 +103,6 @@ void PairTersoffZBL::read_file(char *file)
if (kname == elements[kelement]) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
@ -128,6 +133,11 @@ void PairTersoffZBL::read_file(char *file)
params[nparams].ZBLcut = values.next_double();
params[nparams].ZBLexpscale = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
if (unit_convert) {
params[nparams].biga *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

17
src/MANYBODY/pair_vashishta.cpp
View File

@ -45,6 +45,7 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
nelements = 0;
elements = NULL;
@ -361,9 +362,15 @@ void PairVashishta::read_file(char *file)
// open file on proc 0
if (comm->me == 0) {
PotentialFileReader reader(lmp, file, "Vashishta");
PotentialFileReader reader(lmp, file, "Vashishta", unit_convert_flag);
char * line;
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
while((line = reader.next_line(NPARAMS_PER_LINE))) {
try {
ValueTokenizer values(line);
@ -412,6 +419,14 @@ void PairVashishta::read_file(char *file)
params[nparams].r0 = values.next_double();
params[nparams].bigc = values.next_double();
params[nparams].costheta = values.next_double();
if (unit_convert) {
params[nparams].bigh *= conversion_factor;
params[nparams].bigd *= conversion_factor;
params[nparams].bigw *= conversion_factor;
params[nparams].bigb *= conversion_factor;
}
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
}

14
src/MC/fix_bond_swap.cpp
View File

@ -86,17 +86,11 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp, compute group = fix group
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
std::string cmd = id + std::string("_temp");
id_temp = new char[cmd.size()+1];
strcpy(id_temp,cmd.c_str());
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
modify->add_compute(cmd + " all temp");
tflag = 1;
// initialize atom list

10
src/MC/pair_dsmc.cpp
View File

@ -28,6 +28,7 @@
#include "update.h"
#include "random_mars.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -283,12 +284,9 @@ void PairDSMC::init_style()
celly = (domain->boxhi[1] - domain->boxlo[1])/ncellsy;
cellz = (domain->boxhi[2] - domain->boxlo[2])/ncellsz;
if (comm->me == 0) {
if (screen) fprintf(screen,"DSMC cell size = %g x %g x %g\n",
cellx,celly,cellz);
if (logfile) fprintf(logfile,"DSMC cell size = %g x %g x %g\n",
cellx,celly,cellz);
}
if (comm->me == 0)
utils::logmesg(lmp,fmt::format("DSMC cell size = {} x {} x {}\n",
cellx,celly,cellz));
total_ncells = ncellsx*ncellsy*ncellsz;
vol = cellx*celly*cellz;

10
src/MISC/fix_deposit.cpp
View File

@ -31,6 +31,7 @@
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -286,12 +287,9 @@ void FixDeposit::init()
} else maxradinsert = 0.5;
double separation = MAX(2.0*maxradinsert,maxradall+maxradinsert);
if (sqrt(nearsq) < separation && comm->me == 0) {
char str[128];
sprintf(str,"Fix deposit near setting < possible overlap separation %g",
separation);
error->warning(FLERR,str);
}
if (sqrt(nearsq) < separation && comm->me == 0)
error->warning(FLERR,fmt::format("Fix deposit near setting < possible "
"overlap separation {}",separation));
}
}

23
src/MISC/fix_orient_bcc.cpp
View File

@ -23,6 +23,7 @@
#include <cstring>
#include <cstdlib>
#include <mpi.h>
#include <string>
#include "atom.h"
#include "update.h"
#include "respa.h"
@ -35,6 +36,8 @@
#include "citeme.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -450,20 +453,12 @@ void FixOrientBCC::post_force(int /*vflag*/)
MPI_Allreduce(&maxcount,&max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (screen) fprintf(screen,
"orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
" atoms have %d neighbors\n",
update->ntimestep,atom->natoms,total);
if (logfile) fprintf(logfile,
"orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
" atoms have %d neighbors\n",
update->ntimestep,atom->natoms,total);
if (screen)
fprintf(screen," neighs: min = %d, max = %d, ave = %g\n",
min,max,ave);
if (logfile)
fprintf(logfile," neighs: min = %d, max = %d, ave = %g\n",
min,max,ave);
std::string mesg = fmt::format("orient step {}: {} atoms have {} "
"neighbors\n", update->ntimestep,
atom->natoms,total);
mesg += fmt::format(" neighs: min = {}, max ={}, ave = {}\n",
min,max,ave);
utils::logmesg(lmp,mesg);
}
}
}

23
src/MISC/fix_orient_fcc.cpp
View File

@ -20,6 +20,7 @@
#include <cstring>
#include <cstdlib>
#include <mpi.h>
#include <string>
#include "atom.h"
#include "update.h"
#include "respa.h"
@ -32,6 +33,8 @@
#include "citeme.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -448,20 +451,12 @@ void FixOrientFCC::post_force(int /*vflag*/)
MPI_Allreduce(&maxcount,&max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
if (screen) fprintf(screen,
"orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
" atoms have %d neighbors\n",
update->ntimestep,atom->natoms,total);
if (logfile) fprintf(logfile,
"orient step " BIGINT_FORMAT ": " BIGINT_FORMAT
" atoms have %d neighbors\n",
update->ntimestep,atom->natoms,total);
if (screen)
fprintf(screen," neighs: min = %d, max = %d, ave = %g\n",
min,max,ave);
if (logfile)
fprintf(logfile," neighs: min = %d, max = %d, ave = %g\n",
min,max,ave);
std::string mesg = fmt::format("orient step {}: {} atoms have {} "
"neighbors\n", update->ntimestep,
atom->natoms,total);
mesg += fmt::format(" neighs: min = {}, max ={}, ave = {}\n",
min,max,ave);
utils::logmesg(lmp,mesg);
}
}
}

23
src/MISC/fix_ttm.cpp
View File

@ -29,6 +29,8 @@
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "utils.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -66,25 +68,22 @@ FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
nynodes = force->inumeric(FLERR,arg[11]);
nznodes = force->inumeric(FLERR,arg[12]);
fpr = fopen(arg[13],"r");
if (fpr == NULL) {
char str[128];
snprintf(str,128,"Cannot open file %s",arg[13]);
error->one(FLERR,str);
if (comm->me == 0) {
fpr = fopen(arg[13],"r");
if (fpr == NULL)
error->all(FLERR,fmt::format("Cannot open input file {}: {}",
arg[13], utils::getsyserror()));
}
nfileevery = force->inumeric(FLERR,arg[14]);
if (nfileevery) {
if (narg != 16) error->all(FLERR,"Illegal fix ttm command");
MPI_Comm_rank(world,&me);
if (me == 0) {
if (comm->me == 0) {
fp = fopen(arg[15],"w");
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open fix ttm file %s",arg[15]);
error->one(FLERR,str);
}
if (fp == NULL)
error->one(FLERR,fmt::format("Cannot open output file {}: {}",
arg[15], utils::getsyserror()));
}
}

30
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -36,7 +36,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename):
MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename):
MLIAPDescriptor(lmp)
{
nelements = 0;
@ -154,15 +154,13 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor
void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double **beta, int vflag)
{
int i,j,jnum,ninside;
double delx,dely,delz,evdwl,rsq;
double delx,dely,delz,rsq;
double fij[3];
int *jlist,*numneigh,**firstneigh;
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -245,13 +243,13 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double
// add in global and per-atom virial contributions
// this is optional and has no effect on force calculation
if (vflag)
pairmliap->v_tally(i,j,
fij[0],fij[1],fij[2],
-snaptr->rij[jj][0],-snaptr->rij[jj][1],
-snaptr->rij[jj][2]);
}
}
@ -404,11 +402,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
FILE *fpparam;
if (comm->me == 0) {
fpparam = force->open_potential(paramfilename);
if (fpparam == NULL) {
char str[128];
snprintf(str,128,"Cannot open SNAP parameter file %s",paramfilename);
error->one(FLERR,str);
}
if (fpparam == NULL)
error->one(FLERR,fmt::format("Cannot open SNAP parameter file {}: {}",
paramfilename, utils::getsyserror()));
}
char line[MAXLINE],*ptr;
@ -444,8 +440,8 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
utils::logmesg(lmp, fmt::format("SNAP keyword {} {} \n", keywd, keyval));
}
// check for keywords with one value per element
// check for keywords with one value per element
if (strcmp(keywd,"elems") == 0 ||
strcmp(keywd,"radelems") == 0 ||
strcmp(keywd,"welems") == 0) {
@ -478,9 +474,9 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
} else {
// all other keywords take one value
// all other keywords take one value
if (nwords != 2)
if (nwords != 2)
error->all(FLERR,"Incorrect SNAP parameter file");
if (strcmp(keywd,"nelems") == 0) {
@ -514,8 +510,8 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
}
}
if (!rcutfacflag || !twojmaxflag || !nelementsflag ||
if (!rcutfacflag || !twojmaxflag || !nelementsflag ||
!elementsflag || !radelemflag || !wjelemflag)
error->all(FLERR,"Incorrect SNAP parameter file");

12
src/MLIAP/mliap_model.cpp
View File

@ -22,6 +22,7 @@
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "fmt/format.h"
using namespace LAMMPS_NS;
@ -73,11 +74,9 @@ void MLIAPModel::read_coeffs(char *coefffilename)
FILE *fpcoeff;
if (comm->me == 0) {
fpcoeff = force->open_potential(coefffilename);
if (fpcoeff == NULL) {
char str[128];
snprintf(str,128,"Cannot open MLIAPModel coefficient file %s",coefffilename);
error->one(FLERR,str);
}
if (fpcoeff == NULL)
error->one(FLERR,fmt::format("Cannot open MLIAPModel coeff file {}: {}",
coefffilename,utils::getsyserror()));
}
char line[MAXLINE],*ptr;
@ -153,7 +152,6 @@ void MLIAPModel::read_coeffs(char *coefffilename)
}
if (comm->me == 0) fclose(fpcoeff);
}
/* ----------------------------------------------------------------------
@ -164,9 +162,7 @@ double MLIAPModel::memory_usage()
{
double bytes = 0;
int n = atom->ntypes+1;
bytes += nelements*nparams*sizeof(double); // coeffelem
return bytes;
}

4
src/MLIAP/mliap_model_linear.cpp
View File

@ -28,7 +28,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename)
{
ndescriptors = nparams - 1;
@ -79,7 +79,7 @@ void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double **
for (int icoeff = 0; icoeff < ndescriptors; icoeff++)
etmp += coeffi[icoeff+1]*descriptors[ii][icoeff];
pairmliap->e_tally(i,etmp);
}
}

4
src/MLIAP/mliap_model_quadratic.cpp
View File

@ -28,7 +28,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) :
MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) :
MLIAPModel(lmp, coefffilename)
{
nonlinearflag = 1;
@ -81,7 +81,7 @@ void MLIAPModelQuadratic::gradient(PairMLIAP* pairmliap, NeighList* list, double
// add in contributions to global and per-atom energy
// this is optional and has no effect on force calculation
if (eflag) {
// energy of atom I

4
src/MLIAP/pair_mliap.cpp
View File

@ -93,10 +93,10 @@ void PairMLIAP::compute(int eflag, int vflag)
model->gradient(this, list, descriptors, beta, eflag);
// calculate force contributions beta_i*dB_i/dR_j
descriptor->backward(this, list, beta, vflag);
// calculate stress
// calculate stress
if (vflag_fdotr) virial_fdotr_compute();
}

2
src/MLIAP/pair_mliap.h
View File

@ -45,7 +45,7 @@ protected:
double** beta; // betas for all atoms in list
double** descriptors; // descriptors for all atoms in list
int ndescriptors; // number of descriptors
int ndescriptors; // number of descriptors
int beta_max; // number of atoms allocated for beta, descriptors
class MLIAPModel* model;

Some files were not shown because too many files have changed in this diff Show More