Merge pull request #4181 from gplummer317/ctip
Charge Transfer Interatomic Potential (CTIP) pair style and qeq fix
This commit is contained in:
Axel Kohlmeyer
@ -178,6 +178,7 @@ OPT.
|
||||
* :doc:`python/move <fix_python_move>`
|
||||
* :doc:`qbmsst <fix_qbmsst>`
|
||||
* :doc:`qeq/comb (o) <fix_qeq_comb>`
|
||||
* :doc:`qeq/ctip <fix_qeq>`
|
||||
* :doc:`qeq/dynamic <fix_qeq>`
|
||||
* :doc:`qeq/fire <fix_qeq>`
|
||||
* :doc:`qeq/point <fix_qeq>`
|
||||
|
||||
@ -59,6 +59,7 @@ OPT.
|
||||
* :doc:`comb (o) <pair_comb>`
|
||||
* :doc:`comb3 <pair_comb>`
|
||||
* :doc:`cosine/squared <pair_cosine_squared>`
|
||||
* :doc:`coul/ctip <pair_coul>`
|
||||
* :doc:`coul/cut (gko) <pair_coul>`
|
||||
* :doc:`coul/cut/dielectric <pair_dielectric>`
|
||||
* :doc:`coul/cut/global (o) <pair_coul>`
|
||||
|
||||
@ -357,6 +357,7 @@ accelerated styles exist.
|
||||
* :doc:`python/move <fix_python_move>` - move particles using a Python function during a simulation run
|
||||
* :doc:`qbmsst <fix_qbmsst>` - quantum bath multi-scale shock technique time integrator
|
||||
* :doc:`qeq/comb <fix_qeq_comb>` - charge equilibration for COMB potential
|
||||
* :doc:`qeq/ctip <fix_qeq>` - charge equilibration for CTIP potential
|
||||
* :doc:`qeq/dynamic <fix_qeq>` - charge equilibration via dynamic method
|
||||
* :doc:`qeq/fire <fix_qeq>` - charge equilibration via FIRE minimizer
|
||||
* :doc:`qeq/point <fix_qeq>` - charge equilibration via point method
|
||||
|
||||
@ -1,6 +1,7 @@
|
||||
.. index:: fix qeq/point
|
||||
.. index:: fix qeq/shielded
|
||||
.. index:: fix qeq/slater
|
||||
.. index:: fix qeq/ctip
|
||||
.. index:: fix qeq/dynamic
|
||||
.. index:: fix qeq/fire
|
||||
|
||||
@ -13,6 +14,9 @@ fix qeq/shielded command
|
||||
fix qeq/slater command
|
||||
======================
|
||||
|
||||
fix qeq/ctip command
|
||||
====================
|
||||
|
||||
fix qeq/dynamic command
|
||||
=======================
|
||||
|
||||
@ -27,18 +31,20 @@ Syntax
|
||||
fix ID group-ID style Nevery cutoff tolerance maxiter qfile keyword ...
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* style = *qeq/point* or *qeq/shielded* or *qeq/slater* or *qeq/dynamic* or *qeq/fire*
|
||||
* style = *qeq/point* or *qeq/shielded* or *qeq/slater* or *qeq/ctip* or *qeq/dynamic* or *qeq/fire*
|
||||
* Nevery = perform charge equilibration every this many steps
|
||||
* cutoff = global cutoff for charge-charge interactions (distance unit)
|
||||
* tolerance = precision to which charges will be equilibrated
|
||||
* maxiter = maximum iterations to perform charge equilibration
|
||||
* qfile = a filename with QEq parameters or *coul/streitz* or *reaxff*
|
||||
* qfile = a filename with QEq parameters or *coul/streitz* or *coul/ctip* or *reaxff*
|
||||
* zero or more keyword/value pairs may be appended
|
||||
* keyword = *alpha* or *qdamp* or *qstep* or *warn*
|
||||
* keyword = *alpha* or *cdamp* or *maxrepeat* or *qdamp* or *qstep* or *warn*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*alpha* value = Slater type orbital exponent (qeq/slater only)
|
||||
*cdamp* value = damping parameter for Coulomb interactions (qeq/ctip only)
|
||||
*maxrepeat* value = number of equilibration cycles allowed to ensure no atoms cross charge bounds (qeq/ctip only)
|
||||
*qdamp* value = damping factor for damped dynamics charge solver (qeq/dynamic and qeq/fire only)
|
||||
*qstep* value = time step size for damped dynamics charge solver (qeq/dynamic and qeq/fire only)
|
||||
*warn* value = do (=yes) or do not (=no) print a warning when the maximum number of iterations is reached
|
||||
@ -51,6 +57,7 @@ Examples
|
||||
fix 1 all qeq/point 1 10 1.0e-6 200 param.qeq1
|
||||
fix 1 qeq qeq/shielded 1 8 1.0e-6 100 param.qeq2
|
||||
fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2
|
||||
fix 1 all qeq/ctip 1 12 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
|
||||
fix 1 qeq qeq/dynamic 1 12 1.0e-3 100 my_qeq
|
||||
fix 1 all qeq/fire 1 10 1.0e-3 100 my_qeq qdamp 0.2 qstep 0.1
|
||||
|
||||
@ -103,7 +110,7 @@ equalizes the derivative of energy with respect to charge of all the
|
||||
atoms) by adjusting the partial charge on individual atoms based on
|
||||
interactions with their neighbors within *cutoff*\ . It requires a few
|
||||
parameters in the appropriate units for each atom type which are read
|
||||
from a file specified by *qfile*\ . The file has the following format
|
||||
from a file specified by *qfile*\ . The file has the following format:
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
@ -112,20 +119,32 @@ from a file specified by *qfile*\ . The file has the following format
|
||||
...
|
||||
Ntype chi eta gamma zeta qcore
|
||||
|
||||
except for fix style *qeq/ctip* where the format is:
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
1 chi eta gamma zeta qcore qmin qmax omega
|
||||
2 chi eta gamma zeta qcore qmin qmax omega
|
||||
...
|
||||
Ntype chi eta gamma zeta qcore qmin qmax omega
|
||||
|
||||
There have to be parameters given for every atom type. Wildcard entries
|
||||
are possible using the same type range syntax as for "coeff" commands
|
||||
(i.e., n\*m, n\*, \*m, \*). Later entries will overwrite previous ones.
|
||||
Empty lines or any text following the pound sign (#) are ignored.
|
||||
Each line starts with the atom type followed by five parameters.
|
||||
Only a subset of the parameters is used by each QEq style as described
|
||||
below, thus the others can be set to 0.0 if desired, but all five
|
||||
entries per line are required.
|
||||
Empty lines or any text following the pound sign (#) are ignored. Each
|
||||
line starts with the atom type followed by eight parameters. Only a
|
||||
subset of the parameters is used by each QEq style as described below,
|
||||
thus the others can be set to 0.0 if desired, but all eight entries per
|
||||
line are required.
|
||||
|
||||
* *chi* = electronegativity in energy units
|
||||
* *eta* = self-Coulomb potential in energy units
|
||||
* *gamma* = shielded Coulomb constant defined by :ref:`ReaxFF force field <vanDuin>` in distance units
|
||||
* *zeta* = Slater type orbital exponent defined by the :ref:`Streitz-Mintmire <Streitz1>` potential in reverse distance units
|
||||
* *qcore* = charge of the nucleus defined by the :ref:`Streitz-Mintmire potential <Streitz1>` potential in charge units
|
||||
* *qmin* = lower bound on the allowed charge defined by the :ref:`CTIP <CTIP1>` potential in charge units
|
||||
* *qmax* = upper bound on the allowed charge defined by the :ref:`CTIP <CTIP1>` potential in charge units
|
||||
* *omega* = penalty parameter used to enforce charge bounds defined by the :ref:`CTIP <CTIP1>` potential in energy units
|
||||
|
||||
The fix qeq styles will print a warning if the charges are not
|
||||
equilibrated within *tolerance* by *maxiter* steps, unless the
|
||||
@ -171,6 +190,22 @@ on atoms via the matrix inversion method. A tolerance of 1.0e-6 is
|
||||
usually a good number. Keyword *alpha* can be used to change the Slater
|
||||
type orbital exponent.
|
||||
|
||||
.. versionadded:: TBD
|
||||
|
||||
The *qeq/ctip* style describes partial charges on atoms in the same way
|
||||
as style *qeq/shielded* but also enables the definition of charge
|
||||
bounds. Only the *chi*, *eta*, *gamma*, *qmin*, *qmax*, and *omega*
|
||||
parameters from the *qfile* file are used. When using the string
|
||||
*coul/ctip* as filename, these parameters are extracted directly from an
|
||||
active *coul/ctip* pair style. This style solves partial charges on
|
||||
atoms via the matrix inversion method. Keyword *cdamp* can be used to
|
||||
change the damping parameter used to calculate Coulomb interactions.
|
||||
Keyword *maxrepeat* can be used to adjust the number of equilibration
|
||||
cycles allowed to ensure no atoms have crossed the charge bounds. A
|
||||
value of 10 is usually a good choice. A tolerance between 1.0e-6 and
|
||||
1.0e-8 is usually a good choice but should be checked in conjunction
|
||||
with the timestep for adequate energy conservation during dynamic runs.
|
||||
|
||||
The *qeq/dynamic* style describes partial charges on atoms as point
|
||||
charges that interact through 1/r, but the extended Lagrangian method is
|
||||
used to solve partial charges on atoms. Only the *chi* and *eta*
|
||||
@ -186,7 +221,7 @@ minimization algorithm to solve for equilibrium charges. Keyword
|
||||
*qdamp* can be used to change the damping factor, while keyword *qstep*
|
||||
can be used to change the time step size.
|
||||
|
||||
Note that *qeq/point*, *qeq/shielded*, and *qeq/slater* describe
|
||||
Note that *qeq/point*, *qeq/shielded*, *qeq/slater*, and *qeq/ctip* describe
|
||||
different charge models, whereas the matrix inversion method and the
|
||||
extended Lagrangian method (\ *qeq/dynamic* and *qeq/fire*\ ) are
|
||||
different solvers.
|
||||
@ -266,6 +301,11 @@ Chemistry, 95, 3358-3363 (1991).
|
||||
**(Streitz-Mintmire)** F. H. Streitz, J. W. Mintmire, Physical Review B, 50,
|
||||
16, 11996 (1994)
|
||||
|
||||
.. _CTIP1:
|
||||
|
||||
**(CTIP)** G. Plummer, J. P. Tavenner, M. I. Mendelev, Z. Wu, J. W. Lawson,
|
||||
in preparation
|
||||
|
||||
.. _vanDuin:
|
||||
|
||||
**(ReaxFF)** A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard III, J
|
||||
|
||||
@ -4,6 +4,7 @@
|
||||
.. index:: pair_style coul/cut/omp
|
||||
.. index:: pair_style coul/cut/global
|
||||
.. index:: pair_style coul/cut/global/omp
|
||||
.. index:: pair_style coul/ctip
|
||||
.. index:: pair_style coul/debye
|
||||
.. index:: pair_style coul/debye/gpu
|
||||
.. index:: pair_style coul/debye/kk
|
||||
@ -38,6 +39,9 @@ pair_style coul/cut/global command
|
||||
|
||||
Accelerator Variants: *coul/cut/omp*
|
||||
|
||||
pair_style coul/ctip command
|
||||
============================
|
||||
|
||||
pair_style coul/debye command
|
||||
=============================
|
||||
|
||||
@ -79,7 +83,6 @@ pair_style tip4p/long command
|
||||
|
||||
Accelerator Variants: *tip4p/long/omp*
|
||||
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
@ -87,6 +90,7 @@ Syntax
|
||||
|
||||
pair_style coul/cut cutoff
|
||||
pair_style coul/cut/global cutoff
|
||||
pair_style coul/ctip alpha cutoff
|
||||
pair_style coul/debye kappa cutoff
|
||||
pair_style coul/dsf alpha cutoff
|
||||
pair_style coul/exclude cutoff
|
||||
@ -116,6 +120,9 @@ Examples
|
||||
pair_coeff * *
|
||||
pair_coeff 2 2 3.5
|
||||
|
||||
pair_style coul/ctip 0.30 12.0
|
||||
pair_coeff * * NiO.ctip Ni O
|
||||
|
||||
pair_style coul/debye 1.4 3.0
|
||||
pair_coeff * *
|
||||
pair_coeff 2 2 3.5
|
||||
@ -173,6 +180,32 @@ coulomb styles in :doc:`hybrid pair styles <pair_hybrid>`.
|
||||
|
||||
----------
|
||||
|
||||
.. versionadded:: TBD
|
||||
|
||||
Style *coul/ctip* computes the Coulomb interations as described in
|
||||
:ref:`Plummer <Plummer1>`. It uses the the damped shifted model as in
|
||||
style *coul/dsf* but is further extended to the second derivative of
|
||||
the potential and incorporates empirical charge shielding meant to
|
||||
approximate the more expensive Coulomb integrals used in style *coul/streitz*.
|
||||
More details can be found in the referenced paper. Like the style *coul/streitz*,
|
||||
style *coul/ctip* is a variable charge potential and must be hybridized
|
||||
with a short-range potential via the :doc:`pair_style hybrid/overlay <pair_hybrid>`
|
||||
command. Charge equilibration must be performed with the :doc:`fix qeq/ctip
|
||||
<fix_qeq>` command. For example:
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style hybrid/overlay eam/fs coul/ctip 0.30 12.0
|
||||
pair_coeff * * eam/fs NiO.eam.fs Ni O
|
||||
pair_coeff * * coul/ctip NiO.ctip Ni O
|
||||
fix 1 all qeq/ctip 1 12.0 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
|
||||
|
||||
See the examples/ctip directory for an example input script using the CTIP
|
||||
potential. An Ni-O CTIP and EAM/FS parametrization are included for use with
|
||||
the example.
|
||||
|
||||
----------
|
||||
|
||||
Style *coul/debye* adds an additional exp() damping factor to the
|
||||
Coulombic term, given by
|
||||
|
||||
@ -399,16 +432,18 @@ Restrictions
|
||||
""""""""""""
|
||||
|
||||
The *coul/long*, *coul/msm*, *coul/streitz*, and *tip4p/long* styles are
|
||||
part of the KSPACE package. The *coul/cut/global*, *coul/exclude* styles are
|
||||
part of the EXTRA-PAIR package. The *tip4p/cut* style is part of the MOLECULE
|
||||
package. A pair style is only enabled if LAMMPS was built with its
|
||||
corresponding package. See the :doc:`Build package <Build_package>`
|
||||
doc page for more info.
|
||||
part of the KSPACE package. The *coul/cut/global*, *coul/exclude*, and
|
||||
*coul/ctip* styles are part of the EXTRA-PAIR package. The *tip4p/cut*
|
||||
style is part of the MOLECULE package. A pair style is only enabled if
|
||||
LAMMPS was built with its corresponding package. See the
|
||||
:doc:`Build package <Build_package>` page for more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style, hybrid/overlay <pair_hybrid>`, :doc:`kspace_style <kspace_style>`
|
||||
:doc:`pair_coeff <pair_coeff>`,
|
||||
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
|
||||
:doc:`kspace_style <kspace_style>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
@ -432,6 +467,11 @@ Phys, 110, 8254 (1999).
|
||||
**(Streitz)** F. H. Streitz, J. W. Mintmire, Phys Rev B, 50, 11996-12003
|
||||
(1994).
|
||||
|
||||
.. _Plummer1:
|
||||
|
||||
**(Plummer)** G. Plummer, J. P. Tavenner, M. I. Mendelev, Z. Wu, J. W. Lawson,
|
||||
in preparation
|
||||
|
||||
.. _Jorgensen3:
|
||||
|
||||
**(Jorgensen)** Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem
|
||||
|
||||
@ -151,6 +151,7 @@ accelerated styles exist.
|
||||
* :doc:`comb <pair_comb>` - charge-optimized many-body (COMB) potential
|
||||
* :doc:`comb3 <pair_comb>` - charge-optimized many-body (COMB3) potential
|
||||
* :doc:`cosine/squared <pair_cosine_squared>` - Cooke-Kremer-Deserno membrane model potential
|
||||
* :doc:`coul/ctip <pair_coul>` - Charge Transfer Interatomic (Coulomb) Potential
|
||||
* :doc:`coul/cut <pair_coul>` - cutoff Coulomb potential
|
||||
* :doc:`coul/cut/dielectric <pair_dielectric>` -
|
||||
* :doc:`coul/cut/global <pair_coul>` - cutoff Coulomb potential
|
||||
|
||||
1
examples/streitz/NiO.ctip
Symbolic link
1
examples/streitz/NiO.ctip
Symbolic link
@ -0,0 +1 @@
|
||||
../../potentials/NiO.ctip
|
||||
1
examples/streitz/NiO.eam.fs
Symbolic link
1
examples/streitz/NiO.eam.fs
Symbolic link
@ -0,0 +1 @@
|
||||
../../potentials/NiO.eam.fs
|
||||
1745
examples/streitz/data.ctip
Normal file
1745
examples/streitz/data.ctip
Normal file
File diff suppressed because it is too large
Load Diff
43
examples/streitz/in.ctip
Normal file
43
examples/streitz/in.ctip
Normal file
@ -0,0 +1,43 @@
|
||||
#CTIP potential for NiO
|
||||
#Contributing author: Gabriel Plummer (NASA)
|
||||
|
||||
#Initialize
|
||||
units metal
|
||||
atom_style charge
|
||||
dimension 3
|
||||
boundary p p p
|
||||
|
||||
#Create Structure
|
||||
read_data data.ctip
|
||||
|
||||
#Define Charges
|
||||
group type1 type 1
|
||||
compute charge1 type1 property/atom q
|
||||
compute q1 type1 reduce ave c_charge1
|
||||
group type2 type 2
|
||||
compute charge2 type2 property/atom q
|
||||
compute q2 type2 reduce ave c_charge2
|
||||
|
||||
#Define Potential
|
||||
pair_style hybrid/overlay eam/fs coul/ctip 0.30 12.0
|
||||
pair_coeff * * eam/fs NiO.eam.fs Ni O
|
||||
pair_coeff * * coul/ctip NiO.ctip Ni O
|
||||
fix qeq all qeq/ctip 1 12.0 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
|
||||
|
||||
#Setup
|
||||
timestep 0.001
|
||||
thermo 100
|
||||
thermo_style custom step temp pe lx ly lz pxx pyy pzz c_q1 c_q2
|
||||
|
||||
#Minimization
|
||||
fix relax all box/relax iso 0
|
||||
minimize 1e-10 1e-10 100000 100000
|
||||
unfix relax
|
||||
|
||||
#Dynamics
|
||||
reset_timestep 0
|
||||
variable T equal 1000
|
||||
variable rnd equal round(random(0,999,${T}))
|
||||
velocity all create ${T} ${rnd} mom yes rot yes
|
||||
fix npt all npt temp ${T} ${T} 0.1 iso 0 0 1
|
||||
run 1000
|
||||
186
examples/streitz/log.11Sep24.ctip.g++.1
Normal file
186
examples/streitz/log.11Sep24.ctip.g++.1
Normal file
@ -0,0 +1,186 @@
|
||||
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-269-g5a12c762f3-modified)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
#CTIP potential for NiO
|
||||
#Contributing author: Gabriel Plummer (NASA)
|
||||
|
||||
#Initialize
|
||||
units metal
|
||||
atom_style charge
|
||||
dimension 3
|
||||
boundary p p p
|
||||
|
||||
#Create Structure
|
||||
read_data data.ctip
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (24.719478 24.719478 24.719478)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
1728 atoms
|
||||
read_data CPU = 0.004 seconds
|
||||
|
||||
#Define Charges
|
||||
group type1 type 1
|
||||
864 atoms in group type1
|
||||
compute charge1 type1 property/atom q
|
||||
compute q1 type1 reduce ave c_charge1
|
||||
group type2 type 2
|
||||
864 atoms in group type2
|
||||
compute charge2 type2 property/atom q
|
||||
compute q2 type2 reduce ave c_charge2
|
||||
|
||||
#Define Potential
|
||||
pair_style hybrid/overlay eam/fs coul/ctip 0.30 12.0
|
||||
pair_coeff * * eam/fs NiO.eam.fs Ni O
|
||||
Reading eam/fs potential file NiO.eam.fs with DATE: 2024-04-29
|
||||
pair_coeff * * coul/ctip NiO.ctip Ni O
|
||||
Reading coul/ctip potential file NiO.ctip with DATE: 2024-09-11
|
||||
fix qeq all qeq/ctip 1 12.0 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
|
||||
|
||||
#Setup
|
||||
timestep 0.001
|
||||
thermo 100
|
||||
thermo_style custom step temp pe lx ly lz pxx pyy pzz c_q1 c_q2
|
||||
|
||||
#Minimization
|
||||
fix relax all box/relax iso 0
|
||||
minimize 1e-10 1e-10 100000 100000
|
||||
|
||||
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
|
||||
|
||||
Your simulation uses code contributions which should be cited:
|
||||
|
||||
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
|
||||
|
||||
@Article{Gissinger24,
|
||||
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
|
||||
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
|
||||
journal = {J. Phys. Chem. B},
|
||||
year = 2024,
|
||||
volume = 128,
|
||||
number = 13,
|
||||
pages = {3282–-3297}
|
||||
}
|
||||
|
||||
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
|
||||
|
||||
Neighbor list info ...
|
||||
update: every = 1 steps, delay = 0 steps, check = yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
3 neighbor lists, perpetual/occasional/extra = 3 0 0
|
||||
(1) pair eam/fs, perpetual, trim from (2)
|
||||
attributes: half, newton on, cut 8
|
||||
pair build: trim
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair coul/ctip, perpetual, half/full from (3)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) fix qeq/ctip, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Energy due to 1 extra global DOFs will be included in minimizer energies
|
||||
(src/min.cpp:219)
|
||||
Per MPI rank memory allocation (min/avg/max) = 55.22 | 55.22 | 55.22 Mbytes
|
||||
Step Temp PotEng Lx Ly Lz Pxx Pyy Pzz c_q1 c_q2
|
||||
0 0 -9633.183 24.719478 24.719478 24.719478 -1491.273 -1491.273 -1491.273 1.2374666 -1.2374666
|
||||
6 0 -9633.1929 24.707505 24.707505 24.707505 0.0050470506 0.0050470504 0.0050470502 1.2410908 -1.2410908
|
||||
Loop time of 1.04745 on 1 procs for 6 steps with 1728 atoms
|
||||
|
||||
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = energy tolerance
|
||||
Energy initial, next-to-last, final =
|
||||
-9633.18301850704 -9633.19294329023 -9633.19294333485
|
||||
Force two-norm initial, final = 42.177998 0.00014264694
|
||||
Force max component initial, final = 42.177998 0.00014260857
|
||||
Final line search alpha, max atom move = 0.0079490417 1.1336014e-06
|
||||
Iterations, force evaluations = 6 8
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.3534 | 0.3534 | 0.3534 | 0.0 | 33.74
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.0004205 | 0.0004205 | 0.0004205 | 0.0 | 0.04
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0.69216 | 0.69216 | 0.69216 | 0.0 | 66.08
|
||||
Other | | 0.001461 | | | 0.14
|
||||
|
||||
Nlocal: 1728 ave 1728 max 1728 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 13897 ave 13897 max 13897 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 216000 ave 216000 max 216000 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 2.34317e+06 ave 2.34317e+06 max 2.34317e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2343168
|
||||
Ave neighs/atom = 1356
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
unfix relax
|
||||
|
||||
#Dynamics
|
||||
reset_timestep 0
|
||||
variable T equal 1000
|
||||
variable rnd equal round(random(0,999,${T}))
|
||||
variable rnd equal round(random(0,999,1000))
|
||||
velocity all create ${T} ${rnd} mom yes rot yes
|
||||
velocity all create 1000 ${rnd} mom yes rot yes
|
||||
velocity all create 1000 233 mom yes rot yes
|
||||
fix npt all npt temp ${T} ${T} 0.1 iso 0 0 1
|
||||
fix npt all npt temp 1000 ${T} 0.1 iso 0 0 1
|
||||
fix npt all npt temp 1000 1000 0.1 iso 0 0 1
|
||||
run 1000
|
||||
Per MPI rank memory allocation (min/avg/max) = 54.35 | 54.35 | 54.35 Mbytes
|
||||
Step Temp PotEng Lx Ly Lz Pxx Pyy Pzz c_q1 c_q2
|
||||
0 1000 -9633.1929 24.707505 24.707505 24.707505 15934.991 15754.787 15735.602 1.2410908 -1.2410908
|
||||
100 600.77874 -9528.0602 24.816519 24.816519 24.816519 -16709.313 -15443.072 -17750.832 1.2181988 -1.2181988
|
||||
200 578.84295 -9490.9794 24.812627 24.812627 24.812627 -5879.8769 -4851.4601 -6722.107 1.2254363 -1.2254363
|
||||
300 694.7973 -9478.5512 24.764285 24.764285 24.764285 14639.415 13827.989 13766.766 1.2372201 -1.2372201
|
||||
400 803.93731 -9462.2542 24.866629 24.866629 24.866629 -4644.2854 -6017.2884 -7744.2567 1.2086229 -1.2086229
|
||||
500 893.70492 -9441.0072 24.891756 24.891756 24.891756 -5784.2232 -8219.1644 -4187.392 1.2001224 -1.2001224
|
||||
600 947.21728 -9416.7532 24.863623 24.863623 24.863623 11265.076 12952.469 11331.883 1.2115124 -1.2115124
|
||||
700 1040.4874 -9409.8397 24.933859 24.933859 24.933859 -6570.5927 -8532.7457 -3284.7317 1.1902794 -1.1902794
|
||||
800 1037.9366 -9398.4828 24.935138 24.935138 24.935138 -4681.8676 -2576.4998 -7160.562 1.1880324 -1.1880324
|
||||
900 1049.4211 -9411.115 24.899733 24.899733 24.899733 3250.2695 8121.9271 5945.0542 1.2024086 -1.2024086
|
||||
1000 964.18789 -9412.4405 24.926442 24.926442 24.926442 -7652.603 -5142.9259 -8351.8835 1.1908935 -1.1908935
|
||||
Loop time of 183.932 on 1 procs for 1000 steps with 1728 atoms
|
||||
|
||||
Performance: 0.470 ns/day, 51.092 hours/ns, 5.437 timesteps/s, 9.395 katom-step/s
|
||||
99.7% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 37.347 | 37.347 | 37.347 | 0.0 | 20.30
|
||||
Neigh | 0.22667 | 0.22667 | 0.22667 | 0.0 | 0.12
|
||||
Comm | 0.052514 | 0.052514 | 0.052514 | 0.0 | 0.03
|
||||
Output | 0.00053967 | 0.00053967 | 0.00053967 | 0.0 | 0.00
|
||||
Modify | 146.29 | 146.29 | 146.29 | 0.0 | 79.53
|
||||
Other | | 0.01552 | | | 0.01
|
||||
|
||||
Nlocal: 1728 ave 1728 max 1728 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 13901 ave 13901 max 13901 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 210388 ave 210388 max 210388 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 2.23448e+06 ave 2.23448e+06 max 2.23448e+06 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2234476
|
||||
Ave neighs/atom = 1293.0995
|
||||
Neighbor list builds = 5
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:03:05
|
||||
186
examples/streitz/log.11Sep24.ctip.g++.4
Normal file
186
examples/streitz/log.11Sep24.ctip.g++.4
Normal file
@ -0,0 +1,186 @@
|
||||
LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-269-g5a12c762f3-modified)
|
||||
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
#CTIP potential for NiO
|
||||
#Contributing author: Gabriel Plummer (NASA)
|
||||
|
||||
#Initialize
|
||||
units metal
|
||||
atom_style charge
|
||||
dimension 3
|
||||
boundary p p p
|
||||
|
||||
#Create Structure
|
||||
read_data data.ctip
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (24.719478 24.719478 24.719478)
|
||||
1 by 2 by 2 MPI processor grid
|
||||
reading atoms ...
|
||||
1728 atoms
|
||||
read_data CPU = 0.007 seconds
|
||||
|
||||
#Define Charges
|
||||
group type1 type 1
|
||||
864 atoms in group type1
|
||||
compute charge1 type1 property/atom q
|
||||
compute q1 type1 reduce ave c_charge1
|
||||
group type2 type 2
|
||||
864 atoms in group type2
|
||||
compute charge2 type2 property/atom q
|
||||
compute q2 type2 reduce ave c_charge2
|
||||
|
||||
#Define Potential
|
||||
pair_style hybrid/overlay eam/fs coul/ctip 0.30 12.0
|
||||
pair_coeff * * eam/fs NiO.eam.fs Ni O
|
||||
Reading eam/fs potential file NiO.eam.fs with DATE: 2024-04-29
|
||||
pair_coeff * * coul/ctip NiO.ctip Ni O
|
||||
Reading coul/ctip potential file NiO.ctip with DATE: 2024-09-11
|
||||
fix qeq all qeq/ctip 1 12.0 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10
|
||||
|
||||
#Setup
|
||||
timestep 0.001
|
||||
thermo 100
|
||||
thermo_style custom step temp pe lx ly lz pxx pyy pzz c_q1 c_q2
|
||||
|
||||
#Minimization
|
||||
fix relax all box/relax iso 0
|
||||
minimize 1e-10 1e-10 100000 100000
|
||||
|
||||
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
|
||||
|
||||
Your simulation uses code contributions which should be cited:
|
||||
|
||||
- Type Label Framework: https://doi.org/10.1021/acs.jpcb.3c08419
|
||||
|
||||
@Article{Gissinger24,
|
||||
author = {Jacob R. Gissinger, Ilia Nikiforov, Yaser Afshar, Brendon Waters, Moon-ki Choi, Daniel S. Karls, Alexander Stukowski, Wonpil Im, Hendrik Heinz, Axel Kohlmeyer, and Ellad B. Tadmor},
|
||||
title = {Type Label Framework for Bonded Force Fields in LAMMPS},
|
||||
journal = {J. Phys. Chem. B},
|
||||
year = 2024,
|
||||
volume = 128,
|
||||
number = 13,
|
||||
pages = {3282–-3297}
|
||||
}
|
||||
|
||||
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
|
||||
|
||||
Neighbor list info ...
|
||||
update: every = 1 steps, delay = 0 steps, check = yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 14
|
||||
ghost atom cutoff = 14
|
||||
binsize = 7, bins = 4 4 4
|
||||
3 neighbor lists, perpetual/occasional/extra = 3 0 0
|
||||
(1) pair eam/fs, perpetual, trim from (2)
|
||||
attributes: half, newton on, cut 8
|
||||
pair build: trim
|
||||
stencil: none
|
||||
bin: none
|
||||
(2) pair coul/ctip, perpetual, half/full from (3)
|
||||
attributes: half, newton on
|
||||
pair build: halffull/newton
|
||||
stencil: none
|
||||
bin: none
|
||||
(3) fix qeq/ctip, perpetual
|
||||
attributes: full, newton on
|
||||
pair build: full/bin/atomonly
|
||||
stencil: full/bin/3d
|
||||
bin: standard
|
||||
WARNING: Energy due to 1 extra global DOFs will be included in minimizer energies
|
||||
(src/min.cpp:219)
|
||||
Per MPI rank memory allocation (min/avg/max) = 19.71 | 19.71 | 19.71 Mbytes
|
||||
Step Temp PotEng Lx Ly Lz Pxx Pyy Pzz c_q1 c_q2
|
||||
0 0 -9633.183 24.719478 24.719478 24.719478 -1491.273 -1491.273 -1491.273 1.2374666 -1.2374666
|
||||
6 0 -9633.1929 24.707505 24.707505 24.707505 0.0050498849 0.0050498856 0.0050498855 1.2410908 -1.2410908
|
||||
Loop time of 0.288818 on 4 procs for 6 steps with 1728 atoms
|
||||
|
||||
99.6% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
Minimization stats:
|
||||
Stopping criterion = energy tolerance
|
||||
Energy initial, next-to-last, final =
|
||||
-9633.18301849946 -9633.19294329957 -9633.19294333066
|
||||
Force two-norm initial, final = 42.177998 0.00014286398
|
||||
Force max component initial, final = 42.177998 0.00014268867
|
||||
Final line search alpha, max atom move = 0.0079490643 1.1342414e-06
|
||||
Iterations, force evaluations = 6 8
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.09401 | 0.095886 | 0.097631 | 0.4 | 33.20
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.0024298 | 0.004185 | 0.0060676 | 2.1 | 1.45
|
||||
Output | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Modify | 0.18779 | 0.1878 | 0.18781 | 0.0 | 65.02
|
||||
Other | | 0.000945 | | | 0.33
|
||||
|
||||
Nlocal: 432 ave 432 max 432 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 8593 ave 8593 max 8593 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 54000 ave 54000 max 54000 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
FullNghs: 585792 ave 585792 max 585792 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2343168
|
||||
Ave neighs/atom = 1356
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
unfix relax
|
||||
|
||||
#Dynamics
|
||||
reset_timestep 0
|
||||
variable T equal 1000
|
||||
variable rnd equal round(random(0,999,${T}))
|
||||
variable rnd equal round(random(0,999,1000))
|
||||
velocity all create ${T} ${rnd} mom yes rot yes
|
||||
velocity all create 1000 ${rnd} mom yes rot yes
|
||||
velocity all create 1000 233 mom yes rot yes
|
||||
fix npt all npt temp ${T} ${T} 0.1 iso 0 0 1
|
||||
fix npt all npt temp 1000 ${T} 0.1 iso 0 0 1
|
||||
fix npt all npt temp 1000 1000 0.1 iso 0 0 1
|
||||
run 1000
|
||||
Per MPI rank memory allocation (min/avg/max) = 18.83 | 18.83 | 18.83 Mbytes
|
||||
Step Temp PotEng Lx Ly Lz Pxx Pyy Pzz c_q1 c_q2
|
||||
0 1000 -9633.1929 24.707505 24.707505 24.707505 15934.991 15754.787 15735.602 1.2410908 -1.2410908
|
||||
100 600.77874 -9528.0602 24.816519 24.816519 24.816519 -16709.313 -15443.072 -17750.832 1.2181988 -1.2181988
|
||||
200 578.84295 -9490.9794 24.812627 24.812627 24.812627 -5879.8769 -4851.4601 -6722.107 1.2254363 -1.2254363
|
||||
300 694.7973 -9478.5512 24.764285 24.764285 24.764285 14639.415 13827.989 13766.766 1.2372201 -1.2372201
|
||||
400 803.93731 -9462.2542 24.866629 24.866629 24.866629 -4644.2853 -6017.2884 -7744.2567 1.2086229 -1.2086229
|
||||
500 893.70492 -9441.0072 24.891756 24.891756 24.891756 -5784.2232 -8219.1644 -4187.392 1.2001224 -1.2001224
|
||||
600 947.21728 -9416.7532 24.863623 24.863623 24.863623 11265.076 12952.469 11331.883 1.2115124 -1.2115124
|
||||
700 1040.4874 -9409.8397 24.933859 24.933859 24.933859 -6570.5926 -8532.7456 -3284.7317 1.1902794 -1.1902794
|
||||
800 1037.9366 -9398.4828 24.935138 24.935138 24.935138 -4681.8675 -2576.5001 -7160.5622 1.1880324 -1.1880324
|
||||
900 1049.4211 -9411.115 24.899733 24.899733 24.899733 3250.2695 8121.9274 5945.0541 1.2024086 -1.2024086
|
||||
1000 964.18789 -9412.4405 24.926442 24.926442 24.926442 -7652.603 -5142.926 -8351.8839 1.1908935 -1.1908935
|
||||
Loop time of 63.5378 on 4 procs for 1000 steps with 1728 atoms
|
||||
|
||||
Performance: 1.360 ns/day, 17.649 hours/ns, 15.739 timesteps/s, 27.196 katom-step/s
|
||||
99.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 | 12.737 | 12.782 | 12.839 | 1.1 | 20.12
|
||||
Neigh | 0.082989 | 0.083627 | 0.084898 | 0.3 | 0.13
|
||||
Comm | 0.39773 | 0.45431 | 0.49896 | 5.5 | 0.72
|
||||
Output | 0.00038299 | 0.0004067 | 0.00047658 | 0.0 | 0.00
|
||||
Modify | 50.207 | 50.207 | 50.208 | 0.0 | 79.02
|
||||
Other | | 0.01047 | | | 0.02
|
||||
|
||||
Nlocal: 432 ave 446 max 421 min
|
||||
Histogram: 1 0 1 0 1 0 0 0 0 1
|
||||
Nghost: 8594.5 ave 8608 max 8579 min
|
||||
Histogram: 1 0 0 0 0 1 1 0 0 1
|
||||
Neighs: 52597 ave 54466 max 51209 min
|
||||
Histogram: 1 0 1 0 1 0 0 0 0 1
|
||||
FullNghs: 558619 ave 576961 max 544134 min
|
||||
Histogram: 1 0 1 0 1 0 0 0 0 1
|
||||
|
||||
Total # of neighbors = 2234476
|
||||
Ave neighs/atom = 1293.0995
|
||||
Neighbor list builds = 5
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:01:04
|
||||
4
potentials/NiO.ctip
Normal file
4
potentials/NiO.ctip
Normal file
@ -0,0 +1,4 @@
|
||||
# CTIP potential parameters for Ni-O UNITS: metal DATE: 2024-09-11 CONTRIBUTOR: Gabriel Plummer gabriel.w.plummer@nasa.gov
|
||||
# X (eV) J (eV) gamma (1/ang) zeta (1/ang) Z (e) qmin qmax omega
|
||||
Ni -5.840000 14.157390 0.549290 0.000000 0.000000 0.000000 2.000000 100.000
|
||||
O 0.000000 15.655862 1.428571 0.000000 0.000000 -2.000000 0.000000 100.000
|
||||
18007
potentials/NiO.eam.fs
Normal file
18007
potentials/NiO.eam.fs
Normal file
File diff suppressed because it is too large
Load Diff
@ -90,6 +90,7 @@ cdeam concentration-dependent EAM
|
||||
cgdna potential files for styles in the CG-DNA package
|
||||
comb COMB potential
|
||||
comb3 COMB3 potential
|
||||
ctip Coulombic portion of CTIP potential
|
||||
eam embedded atom method (EAM) single element, DYNAMO funcfl format
|
||||
eam.alloy EAM multi-element alloy, DYNAMO setfl format
|
||||
eam.fs Finnis-Sinclair EAM format (single element or alloy)
|
||||
|
||||
4
src/.gitignore
vendored
4
src/.gitignore
vendored
@ -439,6 +439,8 @@
|
||||
/pair_coul_cut_global.h
|
||||
/pair_coul_streitz.cpp
|
||||
/pair_coul_streitz.h
|
||||
/pair_coul_ctip.cpp
|
||||
/pair_coul_ctip.h
|
||||
/pair_gauss.cpp
|
||||
/pair_gauss.h
|
||||
/pair_lj96_cut.cpp
|
||||
@ -978,6 +980,8 @@
|
||||
/fix_pour.h
|
||||
/fix_qeq_comb.cpp
|
||||
/fix_qeq_comb.h
|
||||
/fix_qeq_ctip.cpp
|
||||
/fix_qeq_ctip.h
|
||||
/fix_qeq_fire.cpp
|
||||
/fix_qeq_fire.h
|
||||
/fix_qeq_reaxff.cpp
|
||||
|
||||
580
src/EXTRA-PAIR/pair_coul_ctip.cpp
Normal file
580
src/EXTRA-PAIR/pair_coul_ctip.cpp
Normal file
@ -0,0 +1,580 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Gabriel Plummer (NASA)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_coul_ctip.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "ewald_const.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "potential_file_reader.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using namespace EwaldConst;
|
||||
|
||||
static constexpr int DELTA = 4;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulCTIP::PairCoulCTIP(LAMMPS *lmp) :
|
||||
Pair(lmp), params(nullptr), shield(nullptr), shieldcu(nullptr), reffc(nullptr),
|
||||
reffcsq(nullptr), reffc4(nullptr), reffc7(nullptr), s2d_shift(nullptr), f_shift(nullptr),
|
||||
e_shift(nullptr), self_factor(nullptr), qeq_x(nullptr), qeq_j(nullptr), qeq_g(nullptr),
|
||||
qeq_z(nullptr), qeq_c(nullptr), qeq_q1(nullptr), qeq_q2(nullptr), qeq_w(nullptr)
|
||||
{
|
||||
single_enable = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulCTIP::~PairCoulCTIP()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
memory->sfree(params);
|
||||
memory->destroy(elem1param);
|
||||
memory->destroy(shield);
|
||||
memory->destroy(shieldcu);
|
||||
memory->destroy(reffc);
|
||||
memory->destroy(reffcsq);
|
||||
memory->destroy(reffc4);
|
||||
memory->destroy(reffc7);
|
||||
memory->destroy(s2d_shift);
|
||||
memory->destroy(f_shift);
|
||||
memory->destroy(e_shift);
|
||||
memory->destroy(self_factor);
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(qeq_x);
|
||||
memory->destroy(qeq_j);
|
||||
memory->destroy(qeq_g);
|
||||
memory->destroy(qeq_z);
|
||||
memory->destroy(qeq_c);
|
||||
memory->destroy(qeq_q1);
|
||||
memory->destroy(qeq_q2);
|
||||
memory->destroy(qeq_w);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::compute(int eflag, int vflag)
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
int iparam_i, jparam_j;
|
||||
double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair;
|
||||
double r, rsq, reff, reffsq, reff4, forcecoul, factor_coul;
|
||||
double prefactor, erfcc, erfcd, t;
|
||||
double selfion;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
selfion = 0.0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double erfcd_cut, t_cut, erfcc_cut;
|
||||
erfcd_cut = exp(-alpha * alpha * cut_coulsq);
|
||||
t_cut = 1.0 / (1.0 + EWALD_P * alpha * cut_coul);
|
||||
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = map[type[i]];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
iparam_i = elem1param[itype];
|
||||
|
||||
selfion = self(¶ms[iparam_i], qtmp);
|
||||
|
||||
if (evflag) ev_tally(i, i, nlocal, 0, 0.0, selfion, 0.0, 0.0, 0.0, 0.0);
|
||||
|
||||
if (eflag) {
|
||||
double e_self = self_factor[iparam_i][iparam_i] * qtmp * qtmp;
|
||||
ev_tally(i, i, nlocal, 0, 0.0, e_self, 0.0, 0.0, 0.0, 0.0);
|
||||
}
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
jtype = map[type[j]];
|
||||
jparam_j = elem1param[jtype];
|
||||
r = sqrt(rsq);
|
||||
reff = cbrt(rsq * r + 1 / shieldcu[iparam_i][jparam_j]);
|
||||
reffsq = reff * reff;
|
||||
reff4 = reffsq * reffsq;
|
||||
prefactor = qqrd2e * qtmp * q[j] / r;
|
||||
erfcd = exp(-alpha * alpha * rsq);
|
||||
t = 1.0 / (1.0 + EWALD_P * alpha * r);
|
||||
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
|
||||
|
||||
forcecoul = prefactor *
|
||||
(erfcc / r + 2.0 * alpha / MY_PIS * erfcd + rsq * r / reff4 - 1 / r -
|
||||
r * f_shift[iparam_i][jparam_j] + r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul)) *
|
||||
r;
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
|
||||
fpair = forcecoul / rsq;
|
||||
|
||||
f[i][0] += delx * fpair;
|
||||
f[i][1] += dely * fpair;
|
||||
f[i][2] += delz * fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
ecoul = prefactor *
|
||||
(erfcc + r / reff - 1 - r * erfcc_cut / cut_coul - r / reffc[iparam_i][jparam_j] +
|
||||
r / cut_coul + r * f_shift[iparam_i][jparam_j] * (r - cut_coul) -
|
||||
r * s2d_shift[iparam_i][jparam_j] * (r - cut_coul) * (r - cut_coul) * 0.5);
|
||||
if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
|
||||
} else
|
||||
ecoul = 0.0;
|
||||
|
||||
if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulCTIP::self(Param *param, double qi)
|
||||
{
|
||||
double s1 = param->chi, s2 = param->eta, s3 = param->qmin, s4 = param->qmax, s5 = param->omega;
|
||||
|
||||
if (qi < s3) {
|
||||
return qi * ((s1 - 2 * s3 * s5) + qi * (0.50 * s2 + s5)) + s3 * s3 * s5;
|
||||
} else if (qi < s4) {
|
||||
return qi * (s1 + qi * (0.50 * s2));
|
||||
} else {
|
||||
return qi * ((s1 - 2 * s4 * s5) + qi * (0.50 * s2 + s5)) + s4 * s4 * s5;
|
||||
}
|
||||
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag, n + 1, n + 1, "pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++) setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
|
||||
|
||||
memory->create(qeq_x, n + 1, "pair:qeq_x");
|
||||
memory->create(qeq_j, n + 1, "pair:qeq_j");
|
||||
memory->create(qeq_g, n + 1, "pair:qeq_g");
|
||||
memory->create(qeq_z, n + 1, "pair:qeq_z");
|
||||
memory->create(qeq_c, n + 1, "pair:qeq_c");
|
||||
memory->create(qeq_q1, n + 1, "pair:qeq_q1");
|
||||
memory->create(qeq_q2, n + 1, "pair:qeq_q2");
|
||||
memory->create(qeq_w, n + 1, "pair:qeq_w");
|
||||
|
||||
map = new int[n + 1];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR, "Illegal pair_style command");
|
||||
|
||||
alpha = utils::numeric(FLERR, arg[0], false, lmp);
|
||||
cut_coul = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::coeff(int narg, char **arg)
|
||||
{
|
||||
if (!allocated) allocate();
|
||||
|
||||
map_element2type(narg - 3, arg + 3);
|
||||
|
||||
// read potential file and initialize potential parameters
|
||||
|
||||
read_file(arg[2]);
|
||||
setup_params();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::init_style()
|
||||
{
|
||||
if (!atom->q_flag) error->all(FLERR, "Pair style coul/ctip requires atom attribute q");
|
||||
|
||||
neighbor->add_request(this);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
|
||||
memory->destroy(shield);
|
||||
memory->destroy(shieldcu);
|
||||
memory->destroy(reffc);
|
||||
memory->destroy(reffcsq);
|
||||
memory->destroy(reffc4);
|
||||
memory->destroy(reffc7);
|
||||
memory->destroy(s2d_shift);
|
||||
memory->destroy(f_shift);
|
||||
memory->destroy(e_shift);
|
||||
memory->destroy(self_factor);
|
||||
|
||||
memory->create(shield, nelements, nelements, "pair:shield");
|
||||
memory->create(shieldcu, nelements, nelements, "pair:shieldcu");
|
||||
memory->create(reffc, nelements, nelements, "pair:reffc");
|
||||
memory->create(reffcsq, nelements, nelements, "pair:reffcsq");
|
||||
memory->create(reffc4, nelements, nelements, "pair:reffc4");
|
||||
memory->create(reffc7, nelements, nelements, "pair:reffc7");
|
||||
memory->create(s2d_shift, nelements, nelements, "pair:s2d_shift");
|
||||
memory->create(f_shift, nelements, nelements, "pair:f_shift");
|
||||
memory->create(e_shift, nelements, nelements, "pair:e_shift");
|
||||
memory->create(self_factor, nelements, nelements, "pair:self_factor");
|
||||
|
||||
double qqrd2e = force->qqrd2e;
|
||||
double cut_coulcu, cut_coul4, alphacu, erfcd_cut, t_cut, erfcc_cut;
|
||||
cut_coulcu = cut_coulsq * cut_coul;
|
||||
cut_coul4 = cut_coulsq * cut_coulsq;
|
||||
alphacu = alpha * alpha * alpha;
|
||||
erfcd_cut = exp(-alpha * alpha * cut_coulsq);
|
||||
t_cut = 1.0 / (1.0 + EWALD_P * alpha * cut_coul);
|
||||
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
|
||||
|
||||
for (int elt1 = 0; elt1 < nelements; elt1++) {
|
||||
for (int elt2 = 0; elt2 < nelements; elt2++) {
|
||||
shield[elt1][elt2] = sqrt(params[elt1].gamma * params[elt2].gamma);
|
||||
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
|
||||
reffc[elt1][elt2] = std::cbrt(cut_coulcu + 1.0 / shieldcu[elt1][elt2]);
|
||||
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
|
||||
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
|
||||
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
|
||||
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coulcu +
|
||||
4.0 * alpha / MY_PIS * erfcd_cut / cut_coulsq + 4.0 * alphacu / MY_PIS * erfcd_cut -
|
||||
2.0 / cut_coulcu + 4.0 * cut_coul4 / reffc7[elt1][elt2] -
|
||||
2.0 * cut_coul / reffc4[elt1][elt2];
|
||||
f_shift[elt1][elt2] = erfcc_cut / cut_coulsq + 2.0 * alpha / MY_PIS * erfcd_cut / cut_coul -
|
||||
1.0 / cut_coulsq + cut_coulsq / reffc4[elt1][elt2];
|
||||
e_shift[elt1][elt2] = 2.0 * erfcc_cut / cut_coul + 2.0 * alpha / MY_PIS * erfcd_cut -
|
||||
2.0 / cut_coul + 1.0 / reffc[elt1][elt2] + cut_coulcu / reffc4[elt1][elt2] +
|
||||
s2d_shift[elt1][elt2] * cut_coulsq * 0.5;
|
||||
self_factor[elt1][elt2] = -(e_shift[elt1][elt2] * 0.5 + alpha / MY_PIS) * qqrd2e;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulCTIP::init_one(int /*i*/, int /*j*/)
|
||||
{
|
||||
return cut_coul;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::read_file(char *file)
|
||||
{
|
||||
memory->sfree(params);
|
||||
params = nullptr;
|
||||
nparams = 0;
|
||||
maxparam = 0;
|
||||
|
||||
// open file on proc 0
|
||||
if (comm->me == 0) {
|
||||
PotentialFileReader reader(lmp, file, "coul/ctip");
|
||||
char *line;
|
||||
|
||||
while ((line = reader.next_line(NPARAMS_PER_LINE))) {
|
||||
try {
|
||||
ValueTokenizer values(line);
|
||||
|
||||
std::string iname = values.next_string();
|
||||
|
||||
// ielement = 1st args
|
||||
int ielement;
|
||||
|
||||
for (ielement = 0; ielement < nelements; ielement++)
|
||||
if (iname == elements[ielement]) break;
|
||||
|
||||
// load up parameter settings and error check their values
|
||||
|
||||
if (nparams == maxparam) {
|
||||
maxparam += DELTA;
|
||||
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
|
||||
|
||||
// make certain all addional allocated storage is initialized
|
||||
// to avoid false positives when checking with valgrind
|
||||
|
||||
memset(params + nparams, 0, DELTA * sizeof(Param));
|
||||
}
|
||||
|
||||
params[nparams].ielement = ielement;
|
||||
params[nparams].chi = values.next_double();
|
||||
params[nparams].eta = values.next_double();
|
||||
params[nparams].gamma = values.next_double();
|
||||
params[nparams].zeta = values.next_double();
|
||||
params[nparams].zcore = values.next_double();
|
||||
params[nparams].qmin = values.next_double();
|
||||
params[nparams].qmax = values.next_double();
|
||||
params[nparams].omega = values.next_double();
|
||||
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
// parameter sanity check
|
||||
|
||||
if (params[nparams].eta < 0.0 || params[nparams].zeta < 0.0 || params[nparams].zcore < 0.0 ||
|
||||
params[nparams].gamma < 0.0 || params[nparams].qmin > params[nparams].qmax ||
|
||||
params[nparams].omega < 0.0)
|
||||
error->one(FLERR, "Illegal coul/ctip parameter");
|
||||
|
||||
nparams++;
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
|
||||
|
||||
if (comm->me != 0) {
|
||||
params = (Param *) memory->srealloc(params, maxparam * sizeof(Param), "pair:params");
|
||||
}
|
||||
|
||||
MPI_Bcast(params, maxparam * sizeof(Param), MPI_BYTE, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::setup_params()
|
||||
{
|
||||
int i, m, n;
|
||||
|
||||
// set elem1param
|
||||
|
||||
memory->destroy(elem1param);
|
||||
memory->create(elem1param, nelements, "pair:elem1param");
|
||||
|
||||
for (i = 0; i < nelements; i++) {
|
||||
n = -1;
|
||||
for (m = 0; m < nparams; m++) {
|
||||
if (i == params[m].ielement) {
|
||||
if (n >= 0) error->all(FLERR, "Potential file has duplicate entry");
|
||||
n = m;
|
||||
}
|
||||
}
|
||||
if (n < 0) error->all(FLERR, "Potential file is missing an entry");
|
||||
elem1param[i] = n;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i, j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j], sizeof(int), 1, fp); }
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i, j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
|
||||
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&alpha, sizeof(double), 1, fp);
|
||||
fwrite(&cut_coul, sizeof(double), 1, fp);
|
||||
fwrite(&offset_flag, sizeof(int), 1, fp);
|
||||
fwrite(&mix_flag, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulCTIP::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR, &alpha, sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&alpha, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairCoulCTIP::extract(const char *str, int &dim)
|
||||
{
|
||||
if (strcmp(str, "cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_coul;
|
||||
}
|
||||
if (strcmp(str, "chi") == 0 && qeq_x) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_x[i] = params[map[i]].chi;
|
||||
else
|
||||
qeq_x[i] = 0.0;
|
||||
return (void *) qeq_x;
|
||||
}
|
||||
if (strcmp(str, "eta") == 0 && qeq_j) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_j[i] = params[map[i]].eta;
|
||||
else
|
||||
qeq_j[i] = 0.0;
|
||||
return (void *) qeq_j;
|
||||
}
|
||||
if (strcmp(str, "gamma") == 0 && qeq_g) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_g[i] = params[map[i]].gamma;
|
||||
else
|
||||
qeq_g[i] = 0.0;
|
||||
return (void *) qeq_g;
|
||||
}
|
||||
if (strcmp(str, "zeta") == 0 && qeq_z) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_z[i] = params[map[i]].zeta;
|
||||
else
|
||||
qeq_z[i] = 0.0;
|
||||
return (void *) qeq_z;
|
||||
}
|
||||
if (strcmp(str, "zcore") == 0 && qeq_c) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_c[i] = params[map[i]].zcore;
|
||||
else
|
||||
qeq_c[i] = 0.0;
|
||||
return (void *) qeq_c;
|
||||
}
|
||||
if (strcmp(str, "qmin") == 0 && qeq_q1) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_q1[i] = params[map[i]].qmin;
|
||||
else
|
||||
qeq_q1[i] = 0.0;
|
||||
return (void *) qeq_q1;
|
||||
}
|
||||
if (strcmp(str, "qmax") == 0 && qeq_q2) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_q2[i] = params[map[i]].qmax;
|
||||
else
|
||||
qeq_q2[i] = 0.0;
|
||||
return (void *) qeq_q2;
|
||||
}
|
||||
if (strcmp(str, "omega") == 0 && qeq_w) {
|
||||
dim = 1;
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
if (map[i] >= 0)
|
||||
qeq_w[i] = params[map[i]].omega;
|
||||
else
|
||||
qeq_w[i] = 0.0;
|
||||
return (void *) qeq_w;
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
70
src/EXTRA-PAIR/pair_coul_ctip.h
Normal file
70
src/EXTRA-PAIR/pair_coul_ctip.h
Normal file
@ -0,0 +1,70 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
// clang-format off
|
||||
PairStyle(coul/ctip,PairCoulCTIP);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_CTIP_H
|
||||
#define LMP_PAIR_COUL_CTIP_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulCTIP : public Pair {
|
||||
public:
|
||||
PairCoulCTIP(class LAMMPS *);
|
||||
~PairCoulCTIP() override;
|
||||
void compute(int, int) override;
|
||||
void settings(int, char **) override;
|
||||
void coeff(int, char **) override;
|
||||
void init_style() override;
|
||||
double init_one(int, int) override;
|
||||
void write_restart(FILE *) override;
|
||||
void read_restart(FILE *) override;
|
||||
void write_restart_settings(FILE *) override;
|
||||
void read_restart_settings(FILE *) override;
|
||||
void *extract(const char *, int &) override;
|
||||
|
||||
static constexpr int NPARAMS_PER_LINE = 9;
|
||||
|
||||
protected:
|
||||
double cut_coul, cut_coulsq;
|
||||
double alpha;
|
||||
|
||||
struct Param {
|
||||
double chi, eta, gamma, zeta, zcore, qmin, qmax, omega;
|
||||
int ielement;
|
||||
};
|
||||
|
||||
Param *params;
|
||||
|
||||
double **shield, **shieldcu, **reffc, **reffcsq, **reffc4, **reffc7;
|
||||
double **s2d_shift, **f_shift, **e_shift, **self_factor;
|
||||
|
||||
double *qeq_x, *qeq_j, *qeq_g, *qeq_z, *qeq_c, *qeq_q1, *qeq_q2, *qeq_w;
|
||||
|
||||
void allocate();
|
||||
|
||||
void read_file(char *);
|
||||
void setup_params();
|
||||
double self(Param *, double);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
@ -44,3 +44,5 @@ action fix_qeq_shielded.cpp
|
||||
action fix_qeq_shielded.h
|
||||
action fix_qeq_slater.cpp
|
||||
action fix_qeq_slater.h
|
||||
action fix_qeq_ctip.cpp
|
||||
action fix_qeq_ctip.h
|
||||
|
||||
@ -27,7 +27,9 @@
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neigh_list.h"
|
||||
#include "pair.h"
|
||||
#include "respa.h"
|
||||
#include "suffix.h"
|
||||
#include "text_file_reader.h"
|
||||
#include "update.h"
|
||||
|
||||
@ -53,7 +55,7 @@ namespace {
|
||||
|
||||
FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), list(nullptr), chi(nullptr), eta(nullptr),
|
||||
gamma(nullptr), zeta(nullptr), zcore(nullptr), chizj(nullptr), shld(nullptr),
|
||||
gamma(nullptr), zeta(nullptr), zcore(nullptr), qmin(nullptr), qmax(nullptr), omega(nullptr), chizj(nullptr), shld(nullptr),
|
||||
s(nullptr), t(nullptr), s_hist(nullptr), t_hist(nullptr), Hdia_inv(nullptr), b_s(nullptr),
|
||||
b_t(nullptr), p(nullptr), q(nullptr), r(nullptr), d(nullptr),
|
||||
qf(nullptr), q1(nullptr), q2(nullptr), qv(nullptr)
|
||||
@ -75,9 +77,9 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
|
||||
error->all(FLERR,"Illegal fix qeq command");
|
||||
|
||||
// must have charges
|
||||
|
||||
if (!atom->q_flag) error->all(FLERR, "Fix {} requires atom attribute q", style);
|
||||
|
||||
alpha = 0.20;
|
||||
swa = 0.0;
|
||||
swb = cutoff;
|
||||
|
||||
@ -115,6 +117,7 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
|
||||
q2 = nullptr;
|
||||
streitz_flag = 0;
|
||||
reax_flag = 0;
|
||||
ctip_flag = 0;
|
||||
qv = nullptr;
|
||||
|
||||
comm_forward = comm_reverse = 1;
|
||||
@ -134,6 +137,8 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
|
||||
streitz_flag = 1;
|
||||
} else if (utils::strmatch(arg[7],"^reax..")) {
|
||||
reax_flag = 1;
|
||||
} else if (strcmp(arg[7],"coul/ctip") == 0) {
|
||||
ctip_flag = 1;
|
||||
} else {
|
||||
read_file(arg[7]);
|
||||
}
|
||||
@ -144,7 +149,6 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixQEq::~FixQEq()
|
||||
{
|
||||
// unregister callbacks to this fix from Atom class
|
||||
|
||||
if (modify->get_fix_by_id(id)) atom->delete_callback(id,Atom::GROW);
|
||||
|
||||
memory->destroy(s_hist);
|
||||
@ -155,12 +159,15 @@ FixQEq::~FixQEq()
|
||||
|
||||
memory->destroy(shld);
|
||||
|
||||
if (!streitz_flag && !reax_flag) {
|
||||
if (!streitz_flag && !reax_flag && !ctip_flag) {
|
||||
memory->destroy(chi);
|
||||
memory->destroy(eta);
|
||||
memory->destroy(gamma);
|
||||
memory->destroy(zeta);
|
||||
memory->destroy(zcore);
|
||||
memory->destroy(qmin);
|
||||
memory->destroy(qmax);
|
||||
memory->destroy(omega);
|
||||
}
|
||||
}
|
||||
|
||||
@ -756,6 +763,9 @@ void FixQEq::read_file(char *file)
|
||||
memory->create(gamma,ntypes+1,"qeq:gamma");
|
||||
memory->create(zeta,ntypes+1,"qeq:zeta");
|
||||
memory->create(zcore,ntypes+1,"qeq:zcore");
|
||||
memory->create(qmin,ntypes+1,"qeq:qmin");
|
||||
memory->create(qmax,ntypes+1,"qeq:qmax");
|
||||
memory->create(omega,ntypes+1,"qeq:omega");
|
||||
|
||||
// read each line out of file, skipping blank lines or leading '#'
|
||||
// store line of params if all 3 element tags are in element list
|
||||
@ -764,7 +774,7 @@ void FixQEq::read_file(char *file)
|
||||
int *setflag = new int[ntypes+1];
|
||||
for (int n=0; n <= ntypes; ++n) {
|
||||
setflag[n] = 0;
|
||||
chi[n] = eta[n] = gamma[n] = zeta[n] = zcore[n] = 0.0;
|
||||
chi[n] = eta[n] = gamma[n] = zeta[n] = zcore[n] = qmin[n] = qmax[n] = omega[n] = 0.0;
|
||||
}
|
||||
|
||||
FILE *fp = nullptr;
|
||||
@ -782,8 +792,13 @@ void FixQEq::read_file(char *file)
|
||||
auto values = reader.next_values(0);
|
||||
|
||||
if (values.count() == 0) continue;
|
||||
if (ctip_flag) {
|
||||
if (values.count() < 9)
|
||||
throw qeq_parser_error(fmt::format("Invalid qeq parameter file for {}", style));
|
||||
} else {
|
||||
if (values.count() < 6)
|
||||
throw qeq_parser_error("Invalid qeq parameter file");
|
||||
throw qeq_parser_error(fmt::format("Invalid qeq parameter file for {}", style));
|
||||
}
|
||||
|
||||
auto word = values.next_string();
|
||||
utils::bounds(FLERR,word,1,ntypes,nlo,nhi,nullptr);
|
||||
@ -800,6 +815,14 @@ void FixQEq::read_file(char *file)
|
||||
for (int n=nlo; n <= nhi; ++n) zeta[n] = val;
|
||||
val = values.next_double();
|
||||
for (int n=nlo; n <= nhi; ++n) zcore[n] = val;
|
||||
if (ctip_flag) {
|
||||
val = values.next_double();
|
||||
for (int n=nlo; n <= nhi; ++n) qmin[n] = val;
|
||||
val = values.next_double();
|
||||
for (int n=nlo; n <= nhi; ++n) qmax[n] = val;
|
||||
val = values.next_double();
|
||||
for (int n=nlo; n <= nhi; ++n) omega[n] = val;
|
||||
}
|
||||
for (int n=nlo; n <= nhi; ++n) setflag[n] = 1;
|
||||
}
|
||||
} catch (EOFException &) {
|
||||
@ -819,4 +842,7 @@ void FixQEq::read_file(char *file)
|
||||
MPI_Bcast(gamma,ntypes+1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(zeta,ntypes+1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(zcore,ntypes+1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(qmin,ntypes+1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(qmax,ntypes+1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(omega,ntypes+1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
|
||||
@ -70,10 +70,10 @@ class FixQEq : public Fix {
|
||||
int maxwarn; // print warning when max iterations was reached
|
||||
double cutoff, cutoff_sq; // neighbor cutoff
|
||||
|
||||
double *chi, *eta, *gamma, *zeta, *zcore; // qeq parameters
|
||||
double *chi, *eta, *gamma, *zeta, *zcore, *qmin, *qmax, *omega;
|
||||
double *chizj;
|
||||
double **shld;
|
||||
int streitz_flag, reax_flag;
|
||||
int streitz_flag, reax_flag, ctip_flag;
|
||||
|
||||
bigint ngroup;
|
||||
|
||||
@ -96,10 +96,6 @@ class FixQEq : public Fix {
|
||||
double *b_s, *b_t;
|
||||
double *p, *q, *r, *d;
|
||||
|
||||
// streitz-mintmire
|
||||
|
||||
double alpha;
|
||||
|
||||
// damped dynamics
|
||||
|
||||
double *qf, *q1, *q2, qdamp, qstep;
|
||||
|
||||
438
src/QEQ/fix_qeq_ctip.cpp
Normal file
438
src/QEQ/fix_qeq_ctip.cpp
Normal file
@ -0,0 +1,438 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Gabriel Plummer (NASA)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_qeq_ctip.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "ewald_const.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "pair.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace EwaldConst;
|
||||
using namespace MathConst;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixQEqCTIP::FixQEqCTIP(LAMMPS *lmp, int narg, char **arg) :
|
||||
FixQEq(lmp, narg, arg), reff(nullptr), reffsq(nullptr), reff4(nullptr), reff7(nullptr),
|
||||
s2d_self(nullptr), shield(nullptr), shieldcu(nullptr), reffc(nullptr), reffcsq(nullptr),
|
||||
reffc4(nullptr), reffc7(nullptr), s2d_shift(nullptr), f_shift(nullptr), e_shift(nullptr)
|
||||
{
|
||||
cdamp = 0.30;
|
||||
maxrepeat = 10;
|
||||
|
||||
// optional args
|
||||
int iarg = 8;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg], "cdamp") == 0) {
|
||||
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip cdamp", error);
|
||||
cdamp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg], "maxrepeat") == 0) {
|
||||
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip maxrepeat", error);
|
||||
maxrepeat = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg], "warn") == 0) {
|
||||
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/ctip warn", error);
|
||||
maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp);
|
||||
iarg += 2;
|
||||
} else
|
||||
error->all(FLERR, "Unknown fix qeq/ctip keyword: {}", arg[iarg]);
|
||||
}
|
||||
|
||||
extract_ctip();
|
||||
}
|
||||
|
||||
FixQEqCTIP::~FixQEqCTIP()
|
||||
{
|
||||
delete[] reff;
|
||||
delete[] reffsq;
|
||||
delete[] reff4;
|
||||
delete[] reff7;
|
||||
delete[] s2d_self;
|
||||
memory->destroy(shield);
|
||||
memory->destroy(shieldcu);
|
||||
memory->destroy(reffc);
|
||||
memory->destroy(reffcsq);
|
||||
memory->destroy(reffc4);
|
||||
memory->destroy(reffc7);
|
||||
memory->destroy(s2d_shift);
|
||||
memory->destroy(f_shift);
|
||||
memory->destroy(e_shift);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::init()
|
||||
{
|
||||
FixQEq::init();
|
||||
|
||||
neighbor->add_request(this, NeighConst::REQ_FULL);
|
||||
|
||||
int ntypes = atom->ntypes;
|
||||
memory->create(shld, ntypes + 1, ntypes + 1, "qeq:shielding");
|
||||
|
||||
delete[] reff;
|
||||
delete[] reffsq;
|
||||
delete[] reff4;
|
||||
delete[] reff7;
|
||||
delete[] s2d_self;
|
||||
|
||||
reff = new double[ntypes];
|
||||
reffsq = new double[ntypes];
|
||||
reff4 = new double[ntypes];
|
||||
reff7 = new double[ntypes];
|
||||
s2d_self = new double[ntypes];
|
||||
|
||||
double r = cutoff;
|
||||
double rsq = r * r;
|
||||
double r6 = rsq * rsq * rsq;
|
||||
|
||||
double erfcd_cut = exp(-cdamp * cdamp * rsq);
|
||||
double t_cut = 1.0 / (1.0 + EWALD_P * cdamp * r);
|
||||
double erfcc_cut =
|
||||
(t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
|
||||
|
||||
for (int elt1 = 0; elt1 < ntypes; elt1++) {
|
||||
reff[elt1] = std::cbrt(rsq * r + 1.0 / (gamma[elt1 + 1] * gamma[elt1 + 1] * gamma[elt1 + 1]));
|
||||
reffsq[elt1] = reff[elt1] * reff[elt1];
|
||||
reff4[elt1] = reffsq[elt1] * reffsq[elt1];
|
||||
reff7[elt1] = reff4[elt1] * reffsq[elt1] * reff[elt1];
|
||||
s2d_self[elt1] = 2.0 * force->qqr2e *
|
||||
(1.5 * erfcc_cut + 2.0 * cdamp / MY_PIS * erfcd_cut +
|
||||
cdamp * cdamp * cdamp / MY_PIS * rsq * erfcd_cut + 0.5 / reff[elt1] - 1.5 / r +
|
||||
r6 / reff7[elt1] + cdamp / MY_PIS);
|
||||
}
|
||||
|
||||
memory->destroy(shield);
|
||||
memory->destroy(shieldcu);
|
||||
memory->destroy(reffc);
|
||||
memory->destroy(reffcsq);
|
||||
memory->destroy(reffc4);
|
||||
memory->destroy(reffc7);
|
||||
memory->destroy(s2d_shift);
|
||||
memory->destroy(f_shift);
|
||||
memory->destroy(e_shift);
|
||||
|
||||
memory->create(shield, ntypes, ntypes, "qeq:shield");
|
||||
memory->create(shieldcu, ntypes, ntypes, "qeq:shieldcu");
|
||||
memory->create(reffc, ntypes, ntypes, "qeq:reffc");
|
||||
memory->create(reffcsq, ntypes, ntypes, "qeq:reffcsq");
|
||||
memory->create(reffc4, ntypes, ntypes, "qeq:reffc4");
|
||||
memory->create(reffc7, ntypes, ntypes, "qeq:reffc7");
|
||||
memory->create(s2d_shift, ntypes, ntypes, "qeq:s2d_shift");
|
||||
memory->create(f_shift, ntypes, ntypes, "qeq:f_shift");
|
||||
memory->create(e_shift, ntypes, ntypes, "qeq:e_shift");
|
||||
|
||||
double cutoffsq = cutoff * cutoff;
|
||||
double cutoffcu = cutoffsq * cutoff;
|
||||
double cutoff4 = cutoffsq * cutoffsq;
|
||||
double cdampcu = cdamp * cdamp * cdamp;
|
||||
|
||||
for (int elt1 = 0; elt1 < atom->ntypes; elt1++) {
|
||||
for (int elt2 = 0; elt2 < atom->ntypes; elt2++) {
|
||||
shield[elt1][elt2] = sqrt(gamma[elt1 + 1] * gamma[elt2 + 1]);
|
||||
shieldcu[elt1][elt2] = shield[elt1][elt2] * shield[elt1][elt2] * shield[elt1][elt2];
|
||||
reffc[elt1][elt2] = cbrt(cutoffcu + 1 / shieldcu[elt1][elt2]);
|
||||
reffcsq[elt1][elt2] = reffc[elt1][elt2] * reffc[elt1][elt2];
|
||||
reffc4[elt1][elt2] = reffcsq[elt1][elt2] * reffcsq[elt1][elt2];
|
||||
reffc7[elt1][elt2] = reffc4[elt1][elt2] * reffcsq[elt1][elt2] * reffc[elt1][elt2];
|
||||
s2d_shift[elt1][elt2] = 2.0 * erfcc_cut / cutoffcu +
|
||||
4.0 * cdamp / MY_PIS * erfcd_cut / cutoffsq + 4.0 * cdampcu / MY_PIS * erfcd_cut -
|
||||
2 / cutoffcu + 4 * cutoff4 / reffc7[elt1][elt2] - 2 * cutoff / reffc4[elt1][elt2];
|
||||
f_shift[elt1][elt2] = erfcc_cut / cutoffsq + 2.0 * cdamp / MY_PIS * erfcd_cut / cutoff -
|
||||
1 / cutoffsq + cutoffsq / reffc4[elt1][elt2];
|
||||
e_shift[elt1][elt2] = erfcc_cut / cutoff + 1 / reffc[elt1][elt2] - 1 / cutoff;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::extract_ctip()
|
||||
{
|
||||
Pair *pair = force->pair_match("^coul/ctip",0);
|
||||
if (pair == nullptr) error->all(FLERR,"No pair style coul/ctip for fix qeq/ctip");
|
||||
int tmp;
|
||||
chi = (double *) pair->extract("chi",tmp);
|
||||
eta = (double *) pair->extract("eta",tmp);
|
||||
gamma = (double *) pair->extract("gamma",tmp);
|
||||
zeta = (double *) pair->extract("zeta",tmp);
|
||||
zcore = (double *) pair->extract("zcore",tmp);
|
||||
qmin = (double *) pair->extract("qmin",tmp);
|
||||
qmax = (double *) pair->extract("qmax",tmp);
|
||||
omega = (double *) pair->extract("omega",tmp);
|
||||
if (chi == nullptr || eta == nullptr || gamma == nullptr || zeta == nullptr ||
|
||||
zcore == nullptr || qmin == nullptr || qmax == nullptr || omega == nullptr)
|
||||
error->all(FLERR, "Fix qeq/ctip could not extract all params from pair style coul/ctip");
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::pre_force(int /*vflag*/)
|
||||
{
|
||||
|
||||
int i,n;
|
||||
|
||||
if (update->ntimestep % nevery) return;
|
||||
|
||||
nlocal = atom->nlocal;
|
||||
|
||||
if (atom->nmax > nmax) reallocate_storage();
|
||||
|
||||
if (nlocal > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE)
|
||||
reallocate_matrix();
|
||||
|
||||
for (i=1; i <= maxrepeat; i++) {
|
||||
init_matvec();
|
||||
matvecs = CG(b_s, s); // CG on s - parallel
|
||||
matvecs += CG(b_t, t); // CG on t - parallel
|
||||
matvecs /= 2;
|
||||
n=calculate_check_Q();
|
||||
MPI_Allreduce(&n, &nout, 1, MPI_INT, MPI_SUM, world);
|
||||
if (nout == 0) break;
|
||||
}
|
||||
|
||||
if (i > maxrepeat && comm->me == 0)
|
||||
error->all(FLERR,"Fix qeq some charges not bound within the domain");
|
||||
|
||||
if (force->kspace) force->kspace->qsum_qsq();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::init_matvec()
|
||||
{
|
||||
compute_H();
|
||||
|
||||
int inum, ii, i;
|
||||
int *ilist;
|
||||
double *q = atom->q, qi;
|
||||
int *type = atom->type;
|
||||
|
||||
double r = cutoff;
|
||||
double rsq = r*r;
|
||||
double r6 = rsq*rsq*rsq;
|
||||
|
||||
double erfcd_cut = exp(-cdamp * cdamp * rsq);
|
||||
double t_cut = 1.0 / (1.0 + EWALD_P * cdamp * r);
|
||||
double erfcc_cut = (t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut) / r;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
|
||||
for (ii = 0; ii < inum; ++ii) {
|
||||
i = ilist[ii];
|
||||
if (atom->mask[i] & groupbit) {
|
||||
|
||||
qi=q[i];
|
||||
if (qi < qmin[type[i]]) {
|
||||
Hdia_inv[i] = 1. / (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1]);
|
||||
b_s[i] = -((chi[type[i]]-2*qmin[type[i]]*omega[type[i]]) + chizj[i]);
|
||||
} else if (qi < qmax[type[i]]) {
|
||||
Hdia_inv[i] = 1. / (eta[type[i]]-s2d_self[type[i]-1]);
|
||||
b_s[i] = -(chi[type[i]] + chizj[i]);
|
||||
} else {
|
||||
Hdia_inv[i] = 1. / (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1]);
|
||||
b_s[i] = -((chi[type[i]]-2*qmax[type[i]]*omega[type[i]]) + chizj[i]);
|
||||
}
|
||||
|
||||
b_t[i] = -1.0;
|
||||
t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]);
|
||||
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
|
||||
}
|
||||
}
|
||||
|
||||
pack_flag = 2;
|
||||
comm->forward_comm(this); //Dist_vector(s);
|
||||
pack_flag = 3;
|
||||
comm->forward_comm(this); //Dist_vector(t);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::compute_H()
|
||||
{
|
||||
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
|
||||
int i, j, ii, jj;
|
||||
double dx, dy, dz, r_sqr, r, reff;
|
||||
double cutoffsq, cutoffcu, cutoff4, cdampcu, erfcc_cut, erfcd_cut, t_cut;
|
||||
double erfcc, erfcd, t;
|
||||
int elt1, elt2;
|
||||
|
||||
double **x = atom->x;
|
||||
int *mask = atom->mask;
|
||||
int *type = atom->type;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
cutoffsq = cutoff * cutoff;
|
||||
cutoffcu = cutoffsq * cutoff;
|
||||
cutoff4 = cutoffsq * cutoffsq;
|
||||
cdampcu = cdamp * cdamp * cdamp;
|
||||
|
||||
erfcd_cut = exp(-cdamp * cdamp * cutoffsq);
|
||||
t_cut = 1.0 / (1.0 + EWALD_P * cdamp * cutoff);
|
||||
erfcc_cut = t_cut * (A1 + t_cut * (A2 + t_cut * (A3 + t_cut * (A4 + t_cut * A5)))) * erfcd_cut;
|
||||
|
||||
|
||||
// fill in the H matrix
|
||||
m_fill = 0;
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
if (mask[i] & groupbit) {
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
H.firstnbr[i] = m_fill;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
dx = x[j][0] - x[i][0];
|
||||
dy = x[j][1] - x[i][1];
|
||||
dz = x[j][2] - x[i][2];
|
||||
r_sqr = dx*dx + dy*dy + dz*dz;
|
||||
|
||||
if (r_sqr <= cutoff_sq) {
|
||||
H.jlist[m_fill] = j;
|
||||
r = sqrt(r_sqr);
|
||||
reff = cbrt(r_sqr * r + 1 / shieldcu[type[i]-1][type[j]-1]);
|
||||
erfcd = exp(-cdamp * cdamp * r_sqr);
|
||||
t = 1.0 / (1.0 + EWALD_P * cdamp * r);
|
||||
erfcc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * erfcd;
|
||||
H.val[m_fill] = 0.5*force->qqr2e*(erfcc/r+1/reff-1/r-e_shift[type[i]-1][type[j]-1]+f_shift[type[i]-1][type[j]-1]*(r-cutoff)-s2d_shift[type[i]-1][type[j]-1]*0.5*(r-cutoff)*(r-cutoff));
|
||||
m_fill++;
|
||||
}
|
||||
}
|
||||
H.numnbrs[i] = m_fill - H.firstnbr[i];
|
||||
}
|
||||
}
|
||||
|
||||
if (m_fill >= H.m)
|
||||
error->all(FLERR,"Fix qeq/ctip has insufficient H matrix size: m_fill={} H.m={}\n",m_fill, H.m);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixQEqCTIP::sparse_matvec(sparse_matrix *A, double *x, double *b)
|
||||
{
|
||||
int i, j, itr_j;
|
||||
double *q=atom->q, qi;
|
||||
int *type = atom->type;
|
||||
|
||||
nlocal = atom->nlocal;
|
||||
nall = atom->nlocal + atom->nghost;
|
||||
|
||||
for (i = 0; i < nlocal; ++i) {
|
||||
if (atom->mask[i] & groupbit) {
|
||||
qi=q[i];
|
||||
if (qi < qmin[type[i]]) {
|
||||
b[i] = (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1])*x[i];
|
||||
} else if (qi < qmax[type[i]]) {
|
||||
b[i] = (eta[type[i]]-s2d_self[type[i]-1]) * x[i];
|
||||
} else {
|
||||
b[i] = (eta[type[i]]+2*omega[type[i]]-s2d_self[type[i]-1])*x[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (i = nlocal; i < nall; ++i) {
|
||||
if (atom->mask[i] & groupbit)
|
||||
b[i] = 0;
|
||||
}
|
||||
|
||||
for (i = 0; i < nlocal; ++i) {
|
||||
if (atom->mask[i] & groupbit) {
|
||||
for (itr_j=A->firstnbr[i]; itr_j<A->firstnbr[i]+A->numnbrs[i]; itr_j++) {
|
||||
j = A->jlist[itr_j];
|
||||
b[i] += A->val[itr_j] * x[j];
|
||||
b[j] += A->val[itr_j] * x[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixQEqCTIP::calculate_check_Q()
|
||||
{
|
||||
int i, k, inum, ii;
|
||||
int *ilist;
|
||||
double u, s_sum, t_sum;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
double qi_old,qi_new;
|
||||
double qi_check1,qi_check2;
|
||||
double qi_check3;
|
||||
int n;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
|
||||
s_sum = parallel_vector_acc( s, inum );
|
||||
t_sum = parallel_vector_acc( t, inum);
|
||||
u = s_sum / t_sum;
|
||||
|
||||
n = 0;
|
||||
for( ii = 0; ii < inum; ++ii ) {
|
||||
i = ilist[ii];
|
||||
if (atom->mask[i] & groupbit) {
|
||||
qi_old = q[i];
|
||||
q[i] = s[i] - u * t[i];
|
||||
|
||||
for( k = 4; k > 0; --k ) {
|
||||
s_hist[i][k] = s_hist[i][k-1];
|
||||
t_hist[i][k] = t_hist[i][k-1];
|
||||
}
|
||||
s_hist[i][0] = s[i];
|
||||
t_hist[i][0] = t[i];
|
||||
|
||||
qi_new = q[i];
|
||||
qi_check1=(qi_new-qmin[type[i]])*(qi_old-qmin[type[i]]);
|
||||
qi_check2=(qi_new-qmax[type[i]])*(qi_old-qmax[type[i]]);
|
||||
if ( qi_check1 < 0.0 || qi_check2 < 0.0 ) {
|
||||
qi_check3=abs(qi_new-qi_old);
|
||||
if (qi_check3 > tolerance) n++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pack_flag = 4;
|
||||
comm->forward_comm( this ); //Dist_vector( atom->q );
|
||||
|
||||
return n;
|
||||
}
|
||||
51
src/QEQ/fix_qeq_ctip.h
Normal file
51
src/QEQ/fix_qeq_ctip.h
Normal file
@ -0,0 +1,51 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
// clang-format off
|
||||
FixStyle(qeq/ctip,FixQEqCTIP);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_QEQ_CTIP_H
|
||||
#define LMP_FIX_QEQ_CTIP_H
|
||||
|
||||
#include "fix_qeq.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixQEqCTIP : public FixQEq {
|
||||
public:
|
||||
FixQEqCTIP(class LAMMPS *, int, char **);
|
||||
~FixQEqCTIP() override;
|
||||
|
||||
void init() override;
|
||||
void pre_force(int) override;
|
||||
|
||||
protected:
|
||||
void init_matvec();
|
||||
void sparse_matvec(sparse_matrix *, double *, double *) override;
|
||||
void compute_H();
|
||||
void extract_ctip();
|
||||
int calculate_check_Q();
|
||||
double *reff, *reffsq, *reff4, *reff7, *s2d_self;
|
||||
double **shield, **shieldcu, **reffc, **reffcsq, **reffc4, **reffc7;
|
||||
double **s2d_shift, **f_shift, **e_shift;
|
||||
|
||||
double cdamp;
|
||||
int maxrepeat;
|
||||
int nout;
|
||||
};
|
||||
} // namespace LAMMPS_NS
|
||||
#endif
|
||||
#endif
|
||||
@ -31,7 +31,7 @@ class FixQEqSlater : public FixQEq {
|
||||
void init() override;
|
||||
void pre_force(int) override;
|
||||
|
||||
private:
|
||||
protected:
|
||||
void init_matvec();
|
||||
void sparse_matvec(sparse_matrix *, double *, double *) override;
|
||||
void compute_H();
|
||||
@ -40,6 +40,7 @@ class FixQEqSlater : public FixQEq {
|
||||
void extract_streitz();
|
||||
|
||||
class PairCoulStreitz *streitz;
|
||||
double alpha;
|
||||
};
|
||||
} // namespace LAMMPS_NS
|
||||
#endif
|
||||
|
||||
92
unittest/force-styles/tests/mol-pair-coul_ctip.yaml
Normal file
92
unittest/force-styles/tests/mol-pair-coul_ctip.yaml
Normal file
@ -0,0 +1,92 @@
|
||||
---
|
||||
lammps_version: 29 Aug 2024
|
||||
tags:
|
||||
date_generated: Tue Sep 10 21:29:25 2024
|
||||
epsilon: 5e-13
|
||||
skip_tests:
|
||||
prerequisites: ! |
|
||||
atom full
|
||||
pair coul/ctip
|
||||
pre_commands: ! |
|
||||
variable units index metal
|
||||
variable newton_pair delete
|
||||
if "$(is_active(package,gpu)) > 0.0" then "variable newton_pair index off" else "variable newton_pair index on"
|
||||
post_commands: ! |
|
||||
timestep 0.001
|
||||
pair_modify mix arithmetic
|
||||
velocity all scale 10.0
|
||||
input_file: in.fourmol
|
||||
pair_style: coul/ctip 0.3 8.0
|
||||
pair_coeff: ! |
|
||||
* * NiO.ctip Ni O Ni O O
|
||||
extract: ! ""
|
||||
natoms: 29
|
||||
init_vdwl: 0
|
||||
init_coul: 268.83312357946403
|
||||
init_stress: ! |2-
|
||||
6.3490063021156979e+00 3.7408695241730894e+00 7.0361156620229268e+00 2.0792271323016851e+00 1.5573754855973017e+00 2.7470162844995999e-01
|
||||
init_forces: ! |2
|
||||
1 -6.6766498869476165e-01 -5.5042307531553758e-01 9.7032502239632901e-01
|
||||
2 7.7521977879752624e-01 5.0872136823803349e-01 -9.8588006812149209e-01
|
||||
3 7.0494670150651594e-03 6.7222063121221878e-02 2.9642961180318752e-02
|
||||
4 -2.1429412404391085e-02 -2.0950579072655018e-02 -3.3568909100718497e-02
|
||||
5 -3.7363529861740602e-02 -2.4923337785193594e-02 1.1516900374344946e-02
|
||||
6 -1.0946612479959126e+00 8.9700509334645628e-01 4.9320185843336117e-01
|
||||
7 1.1854416475916033e-01 -2.6576678610813409e-01 -1.0899734140877935e+00
|
||||
8 6.4955585415778372e-01 -7.6274120860081573e-01 -6.3166430197088130e-01
|
||||
9 3.3316295347683439e-01 8.0330773793723831e-02 1.3555612385253035e+00
|
||||
10 1.1467966540573368e-01 -2.8828389443679403e-01 3.8986752595171752e-02
|
||||
11 -2.4384183205874675e-02 3.3909937515907505e-02 1.8524465479279953e-02
|
||||
12 -1.9578368735092340e-02 -6.4486478848736154e-04 1.4237271164103524e-01
|
||||
13 1.4680378874991862e-01 -6.6502911869983319e-02 -6.2810542877883603e-03
|
||||
14 -1.0152770779985132e-01 2.4657561691070650e-02 -1.6160131249611986e-01
|
||||
15 2.0317060148768745e-02 1.5135710002923145e-01 8.7520151900347225e-03
|
||||
16 4.0894956691719275e-01 -3.8805473224896847e-01 -1.1465067428171507e+00
|
||||
17 -5.6419218424954654e-01 5.9312634326499702e-01 9.1123225799149365e-01
|
||||
18 -1.0353302404272195e-01 -3.0606707755299051e-01 1.4428249631391288e+00
|
||||
19 -7.4737534024659158e-01 -5.4784645392478748e-01 -8.3742536745705298e-01
|
||||
20 8.2767150614670815e-01 8.5638074072256409e-01 -5.2830130487714022e-01
|
||||
21 -4.2535085192871286e-01 -4.2523240078604435e-01 1.4034119287909781e+00
|
||||
22 -7.0318546735293863e-01 -1.9656967920669050e-01 -1.0670964282227489e+00
|
||||
23 1.1119704560782739e+00 6.3955433504897763e-01 -3.1639917324870870e-01
|
||||
24 2.9868496682933654e-01 -1.2313751450902117e+00 7.0609939414552292e-01
|
||||
25 -1.0090099351955220e+00 1.1559130801748081e-01 -8.2798310388092733e-01
|
||||
26 6.9019280914785663e-01 1.1080935733374833e+00 1.0138627381567281e-01
|
||||
27 1.9301178235023941e-01 -1.3471486311081766e+00 4.8386252002863672e-01
|
||||
28 -1.0539921061236890e+00 4.0997619288717529e-01 -7.0289214632426211e-01
|
||||
29 8.7743452785694820e-01 9.3660438688114711e-01 2.1787206316617308e-01
|
||||
run_vdwl: 0
|
||||
run_coul: 268.4476086884773
|
||||
run_stress: ! |2-
|
||||
5.9299248883876601e+00 3.4781130547249339e+00 6.5801272741781647e+00 1.8619728227105692e+00 1.4659099790426291e+00 3.5924356392720608e-01
|
||||
run_forces: ! |2
|
||||
1 -6.0427571725450790e-01 -4.9892889777278226e-01 8.8545319536786815e-01
|
||||
2 7.1312454132801029e-01 4.5923931150828301e-01 -9.0158273806646261e-01
|
||||
3 6.5930657976885961e-03 6.8086911518340634e-02 3.0298419004164478e-02
|
||||
4 -2.0005152917110077e-02 -2.0495734577847459e-02 -3.3168698480578490e-02
|
||||
5 -3.6989504598576979e-02 -2.4815949113475184e-02 1.0615514262009192e-02
|
||||
6 -1.0766475299916556e+00 8.8359720904052264e-01 4.7219323578273070e-01
|
||||
7 1.2279399227984265e-01 -2.6624889781330596e-01 -1.0434209735871274e+00
|
||||
8 6.5311135505192497e-01 -7.1476474929937395e-01 -5.5684953652533686e-01
|
||||
9 3.1104308514004397e-01 3.7965971349324459e-02 1.2513806890774151e+00
|
||||
10 1.1744727580125273e-01 -2.9036328040419229e-01 3.9026787957685875e-02
|
||||
11 -2.4486616434667621e-02 3.4101196287017937e-02 1.7731496391806937e-02
|
||||
12 -1.8112753214851833e-02 -1.5329987398169501e-03 1.3724246989467195e-01
|
||||
13 1.4421485545442478e-01 -6.4023023867411502e-02 -6.4830241578088810e-03
|
||||
14 -1.0093704405052448e-01 2.3626240487656585e-02 -1.5594488783666602e-01
|
||||
15 1.8113350861168019e-02 1.5003134533333243e-01 9.0410683517616990e-03
|
||||
16 3.9564305615954043e-01 -3.8326160090034744e-01 -1.1107744422444059e+00
|
||||
17 -5.5350195813896730e-01 5.9523443686042754e-01 8.7501993236997910e-01
|
||||
18 -3.5281227598702475e-02 -1.9465990896490087e-01 1.2667234231080493e+00
|
||||
19 -6.8038617394630396e-01 -5.1379378261297537e-01 -7.5333352330996883e-01
|
||||
20 6.8989560916162340e-01 7.1158296007868249e-01 -4.3113760667659928e-01
|
||||
21 -3.8916912422204580e-01 -3.6027897454313712e-01 1.2562100641596712e+00
|
||||
22 -6.3524203874677565e-01 -1.8438227441572771e-01 -9.5303869389459417e-01
|
||||
23 1.0072567206938428e+00 5.6287061853542042e-01 -2.8243135805657577e-01
|
||||
24 2.9449971576480227e-01 -1.1247166291512272e+00 6.5352065873792231e-01
|
||||
25 -9.4482316277239342e-01 9.9894478313402951e-02 -7.7751414482791736e-01
|
||||
26 6.2907073683814840e-01 1.0167291975732295e+00 1.0258032074544736e-01
|
||||
27 1.7910891405692830e-01 -1.2099165608953455e+00 4.2142835781775040e-01
|
||||
28 -9.4665562310478557e-01 3.6784596735523806e-01 -6.2528460074895109e-01
|
||||
29 7.8459735260262686e-01 8.4137741883098838e-01 2.0249859538405895e-01
|
||||
...
|
||||
Reference in New Issue
Block a user