diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 55445bfc33..04d1a9969a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -187,6 +187,7 @@ OPT. * :doc:`qeq/slater ` * :doc:`qmmm ` * :doc:`qtb ` + * :doc:`qtpie/reaxff ` * :doc:`rattle ` * :doc:`reaxff/bonds (k) ` * :doc:`reaxff/species (k) ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index ee52be224e..9af607601b 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -366,6 +366,7 @@ accelerated styles exist. * :doc:`qeq/slater ` - charge equilibration via Slater method * :doc:`qmmm ` - functionality to enable a quantum mechanics/molecular mechanics coupling * :doc:`qtb ` - implement quantum thermal bath scheme +* :doc:`qtpie/reaxff ` - apply QTPIE charge equilibration * :doc:`rattle ` - RATTLE constraints on bonds and/or angles * :doc:`reaxff/bonds ` - write out ReaxFF bond information * :doc:`reaxff/species ` - write out ReaxFF molecule information diff --git a/doc/src/fix_acks2_reaxff.rst b/doc/src/fix_acks2_reaxff.rst index ebb1b48051..79a9cf8ea6 100644 --- a/doc/src/fix_acks2_reaxff.rst +++ b/doc/src/fix_acks2_reaxff.rst @@ -111,7 +111,8 @@ LAMMPS was built with that package. See the :doc:`Build package This fix does not correctly handle interactions involving multiple periodic images of the same atom. Hence, it should not be used for -periodic cell dimensions less than :math:`10~\AA`. +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, @@ -122,7 +123,8 @@ components in non-periodic directions. Related commands """""""""""""""" -:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff ` +:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff `, +:doc:`fix qtpi/reaxff ` Default """"""" diff --git a/doc/src/fix_qeq_reaxff.rst b/doc/src/fix_qeq_reaxff.rst index e90842ea6a..e1a09c4fc3 100644 --- a/doc/src/fix_qeq_reaxff.rst +++ b/doc/src/fix_qeq_reaxff.rst @@ -124,7 +124,8 @@ LAMMPS was built with that package. See the :doc:`Build package This fix does not correctly handle interactions involving multiple periodic images of the same atom. Hence, it should not be used for -periodic cell dimensions less than 10 Angstroms. +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. This fix may be used in combination with :doc:`fix efield ` and will apply the external electric field during charge equilibration, @@ -138,7 +139,8 @@ as an atom-style variable using the *potential* keyword for `fix efield`. Related commands """""""""""""""" -:doc:`pair_style reaxff `, :doc:`fix qeq/shielded ` +:doc:`pair_style reaxff `, :doc:`fix qeq/shielded `, +:doc:`fix acks2/reaxff `, :doc:`fix qtpie/reaxff ` Default """"""" diff --git a/doc/src/fix_qtpie_reaxff.rst b/doc/src/fix_qtpie_reaxff.rst new file mode 100644 index 0000000000..cf59e4701a --- /dev/null +++ b/doc/src/fix_qtpie_reaxff.rst @@ -0,0 +1,198 @@ +.. index:: fix qtpie/reaxff + +fix qtpie/reaxff command +======================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + fix ID group-ID qtpie/reaxff Nevery cutlo cuthi tolerance params gfile args + +* ID, group-ID are documented in :doc:`fix ` command +* qtpie/reaxff = style name of this fix command +* Nevery = perform QTPIE every this many steps +* cutlo,cuthi = lo and hi cutoff for Taper radius +* tolerance = precision to which charges will be equilibrated +* params = reaxff or a filename +* gfile = the name of a file containing Gaussian orbital exponents +* one or more keywords or keyword/value pairs may be appended + + .. parsed-literal:: + + keyword = *maxiter* + *maxiter* N = limit the number of iterations to *N* + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff exp.qtpie + fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 params.qtpie exp.qtpie maxiter 500 + +Description +""""""""""" + +The QTPIE charge equilibration method is an extension of the QEq charge +equilibration method. With QTPIE, the partial charges on individual atoms +are computed by minimizing the electrostatic energy of the system in the +same way as the QEq method but where the absolute electronegativity, +:math:`\chi_i`, of each atom in the QEq charge equilibration scheme +:ref:`(Rappe and Goddard) ` is replaced with an effective +electronegativity given by :ref:`(Chen) ` + +.. math:: + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j) S_{ij}} + {\sum_{m=1}^{N}S_{im}}, + +which acts to penalize long-range charge transfer seen with the QEq charge +equilibration scheme. In this equation, :math:`N` is the number of atoms in +the system and :math:`S_{ij}` is the overlap integral between atom :math:`i` +and atom :math:`j`. + +The effect of an external electric field can be incorporated into the QTPIE +method by modifying the absolute or effective electronegativities of each +atom :ref:`(Chen) `. This fix models the effect of an external +electric field by using the effective electronegativity given in +:ref:`(Gergs) `: + +.. math:: + \chi_{\mathrm{eff},i} = \frac{\sum_{j=1}^{N} (\chi_i - \chi_j + \phi_i - \phi_j) S_{ij}} + {\sum_{m=1}^{N}S_{im}}, + +where :math:`\phi_i` and :math:`\phi_j` are the electric +potentials at the positions of atom :math:`i` and :math:`j` +due to the external electric field. + +This fix is typically used in conjunction with the ReaxFF force +field model as implemented in the :doc:`pair_style reaxff ` +command, but it can be used with any potential in LAMMPS, so long as it +defines and uses charges on each atom. For more technical details about the +charge equilibration performed by `fix qtpie/reaxff`, which is the same as in +:doc:`fix qeq/reaxff ` except for the use of +:math:`\chi_{\mathrm{eff},i}`, please refer to :ref:`(Aktulga) `. +To be explicit, this fix replaces :math:`\chi_k` of eq. 3 in +:ref:`(Aktulga) ` with :math:`\chi_{\mathrm{eff},k}`. + +This fix requires the absolute electronegativity, :math:`\chi`, in eV, the +self-Coulomb potential, :math:`\eta`, in eV, and the shielded Coulomb +constant, :math:`\gamma`, in :math:`\AA^{-1}`. If the *params* setting above +is the word "reaxff", then these are extracted from the +:doc:`pair_style reaxff ` command and the ReaxFF force field +file it reads in. If a file name is specified for *params*, then the +parameters are taken from the specified file and the file must contain +one line for each atom type. The latter form must be used when performing +QTPIE with a non-ReaxFF potential. Each line should be formatted as follows, +ensuring that the parameters are given in units of eV, eV, and :math:`\AA^{-1}`, +respectively: + +.. parsed-literal:: + + itype chi eta gamma + +where *itype* is the atom type from 1 to Ntypes. Note that eta is +defined here as twice the eta value in the ReaxFF file. + +The overlap integrals in the equation for :math:`\chi_{\mathrm{eff},i}` +are computed by using normalized 1s Gaussian type orbitals. The Gaussian +orbital exponents, :math:`\alpha`, that are needed to compute the overlap +integrals are taken from the file given by *gfile*. +This file must contain one line for each atom type and provide the Gaussian +orbital exponent for each atom type in units of inverse square Bohr radius. +Each line should be formatted as follows: + +.. parsed-literal:: + + itype alpha + +Empty lines or any text following the pound sign (#) are ignored. An example +*gfile* for a system with two atom types is + +.. parsed-literal:: + + # An example gfile. Exponents are taken from Table 2.2 of Chen, J. (2009). + # Theory and applications of fluctuating-charge models. + # The units of the exponents are 1 / (Bohr radius)^2 . + 1 0.2240 # O + 2 0.5434 # H + +The optional *maxiter* keyword allows changing the max number +of iterations in the linear solver. The default value is 200. + +.. note:: + + In order to solve the self-consistent equations for electronegativity + equalization, LAMMPS imposes the additional constraint that all the + charges in the fix group must add up to zero. The initial charge + assignments should also satisfy this constraint. LAMMPS will print a + warning if that is not the case. + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files +`. This fix computes a global scalar (the number of +iterations) and a per-atom vector (the effective electronegativity), which +can be accessed by various :doc:`output commands `. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. + +This fix is invoked during :doc:`energy minimization `. + +Restrictions +"""""""""""" + +This fix is part of the REAXFF package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +This fix does not correctly handle interactions involving multiple +periodic images of the same atom. Hence, it should not be used for +periodic cell dimensions smaller than the non-bonded cutoff radius, +which is typically :math:`10~\AA` for ReaxFF simulations. + +This fix may be used in combination with :doc:`fix efield ` +and will apply the external electric field during charge equilibration, +but there may be only one fix efield instance used and the electric field +must be applied to all atoms in the system. Consequently, `fix efield` must +be used with *group-ID* all and must not be used with the keyword *region*. +Equal-style variables can be used for electric field vector +components without any further settings. Atom-style variables can be used +for spatially-varying electric field vector components, but the resulting +electric potential must be specified as an atom-style variable using +the *potential* keyword for `fix efield`. + +Related commands +"""""""""""""""" + +:doc:`pair_style reaxff `, :doc:`fix qeq/reaxff `, +:doc:`fix acks2/reaxff ` + +Default +""""""" + +maxiter 200 + +---------- + +.. _Rappe3: + +**(Rappe)** Rappe and Goddard III, Journal of Physical Chemistry, 95, +3358-3363 (1991). + +.. _qtpie-Chen: + +**(Chen)** Chen, Jiahao. Theory and applications of fluctuating-charge models. +University of Illinois at Urbana-Champaign, 2009. + +.. _Gergs: + +**(Gergs)** Gergs, Dirkmann and Mussenbrock. +Journal of Applied Physics 123.24 (2018). + +.. _qeq-Aktulga2: + +**(Aktulga)** Aktulga, Fogarty, Pandit, Grama, Parallel Computing, 38, +245-259 (2012). diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index 03d53d1ff4..495572dc0e 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -20,7 +20,7 @@ Syntax .. parsed-literal:: keyword = *checkqeq* or *lgvdw* or *safezone* or *mincap* or *minhbonds* or *tabulate* or *list/blocking* - *checkqeq* value = *yes* or *no* = whether or not to require qeq/reaxff or acks2/reaxff fix + *checkqeq* value = *yes* or *no* = whether or not to require one of fix qeq/reaxff, fix acks2/reaxff or fix qtpie/reaxff *enobonds* value = *yes* or *no* = whether or not to tally energy of atoms with no bonds *lgvdw* value = *yes* or *no* = whether or not to use a low gradient vdW correction *safezone* = factor used for array allocation @@ -120,20 +120,22 @@ up that process. The ReaxFF parameter files provided were created using a charge equilibration (QEq) model for handling the electrostatic interactions. -Therefore, by default, LAMMPS requires that either the -:doc:`fix qeq/reaxff ` or the -:doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` -command be used with -*pair_style reaxff* when simulating a ReaxFF model, to equilibrate -the charges each timestep. +Therefore, by default, LAMMPS requires that +:doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` +or :doc:`fix acks2/reaxff ` +or :doc:`fix qtpie/reaxff ` +is used with *pair_style reaxff* when simulating a ReaxFF model, +to equilibrate the charges at each timestep. +See the :doc:`fix qeq/reaxff ` or :doc:`fix qeq/shielded ` +or :doc:`fix acks2/reaxff ` +or :doc:`fix qtpie/reaxff ` +command documentation for more details. Using the keyword *checkqeq* with the value *no* turns off the check for the QEq fixes, allowing a simulation to be run without charge equilibration. In this case, the static charges you assign to each atom will be used for computing the electrostatic interactions in -the system. See the :doc:`fix qeq/reaxff ` or -:doc:`fix qeq/shielded ` or :doc:`fix acks2/reaxff ` -command documentation for more details. +the system. Using the optional keyword *lgvdw* with the value *yes* turns on the low-gradient correction of ReaxFF for long-range London Dispersion, @@ -372,8 +374,8 @@ Related commands """""""""""""""" :doc:`pair_coeff `, :doc:`fix qeq/reaxff `, -:doc:`fix acks2/reaxff `, :doc:`fix reaxff/bonds `, -:doc:`fix reaxff/species `, +:doc:`fix acks2/reaxff `, :doc:`fix qtpie/reaxff `, +:doc:`fix reaxff/bonds `, :doc:`fix reaxff/species `, :doc:`compute reaxff/atom ` Default diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 0708a14882..8e601d6c16 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -815,6 +815,7 @@ dipoleflag dir Direc dirname +Dirkmann discoverable discretization discretized @@ -976,6 +977,7 @@ elaplong elastance Electroneg electronegative +electronegativities electronegativity electroneutral electroneutrality @@ -1293,6 +1295,7 @@ Geocomputing georg Georg Geotechnica +Gergs germain Germann Germano @@ -1305,6 +1308,7 @@ gettimeofday geturl gewald Gezelter +gfile Gflop gfortran ghostneigh @@ -1712,6 +1716,7 @@ Jewett jgissing ji Jiang +Jiahao Jiao jik JIK @@ -2366,6 +2371,7 @@ mui Mukherjee Mulders Müller +Mulliken mult multi multibody @@ -2390,6 +2396,7 @@ Murdick Murtola Murty Muser +Mussenbrock mutexes Muto muVT @@ -3073,6 +3080,7 @@ qqr qqrd Qsb qtb +qtpie quadratically quadrupolar quadrupole diff --git a/examples/reaxff/water/gauss_exp.txt b/examples/reaxff/water/gauss_exp.txt new file mode 100644 index 0000000000..4210471e9f --- /dev/null +++ b/examples/reaxff/water/gauss_exp.txt @@ -0,0 +1,5 @@ +# Gaussian orbital exponents (required for fix qtpie/reaxff) taken from Table 2.2 +# of Chen, J. (2009). Theory and applications of fluctuating-charge models. +# The units of the exponents are 1 / (Bohr radius)^2 . +1 0.2240 # O +2 0.5434 # H diff --git a/examples/reaxff/water/in.water.qtpie b/examples/reaxff/water/in.water.qtpie new file mode 100644 index 0000000000..a8f8759444 --- /dev/null +++ b/examples/reaxff/water/in.water.qtpie @@ -0,0 +1,29 @@ +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 diff --git a/examples/reaxff/water/in.water.qtpie.field b/examples/reaxff/water/in.water.qtpie.field new file mode 100644 index 0000000000..e5ac77484f --- /dev/null +++ b/examples/reaxff/water/in.water.qtpie.field @@ -0,0 +1,30 @@ +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 new file mode 100644 index 0000000000..33221ff080 --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.1 @@ -0,0 +1,127 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.056 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + 3000 atoms + replicate CPU = 0.001 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +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} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +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 = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes + Step Temp Press Density Volume + 0 300 10137.041 1 29915.273 + 10 296.09128 3564.7969 1 29915.273 + 20 293.04308 10299.201 1 29915.273 +Loop time of 10.7863 on 1 procs for 20 steps with 3000 atoms + +Performance: 0.080 ns/day, 299.620 hours/ns, 1.854 timesteps/s, 5.563 katom-step/s +100.0% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7275 | 4.7275 | 4.7275 | 0.0 | 43.83 +Neigh | 0.17533 | 0.17533 | 0.17533 | 0.0 | 1.63 +Comm | 0.0017376 | 0.0017376 | 0.0017376 | 0.0 | 0.02 +Output | 8.2065e-05 | 8.2065e-05 | 8.2065e-05 | 0.0 | 0.00 +Modify | 5.8812 | 5.8812 | 5.8812 | 0.0 | 54.52 +Other | | 0.0005226 | | | 0.00 + +Nlocal: 3000 ave 3000 max 3000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 11077 ave 11077 max 11077 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 971775 ave 971775 max 971775 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 971775 +Ave neighs/atom = 323.925 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:12 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 new file mode 100644 index 0000000000..07a348604e --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie-field.g++.4 @@ -0,0 +1,127 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.053 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + 3000 atoms + replicate CPU = 0.002 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 +fix 3 all efield 0.0 0.0 0.05 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +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} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +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 = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes + Step Temp Press Density Volume + 0 300 10137.041 1 29915.273 + 10 296.09128 3564.7969 1 29915.273 + 20 293.04308 10299.201 1 29915.273 +Loop time of 3.14492 on 4 procs for 20 steps with 3000 atoms + +Performance: 0.275 ns/day, 87.359 hours/ns, 6.359 timesteps/s, 19.078 katom-step/s +99.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 | 1.6557 | 1.6847 | 1.7281 | 2.1 | 53.57 +Neigh | 0.086503 | 0.086968 | 0.087627 | 0.2 | 2.77 +Comm | 0.003309 | 0.046699 | 0.075729 | 12.4 | 1.48 +Output | 5.0156e-05 | 5.483e-05 | 6.8111e-05 | 0.0 | 0.00 +Modify | 1.3254 | 1.3261 | 1.3266 | 0.0 | 42.16 +Other | | 0.0004552 | | | 0.01 + +Nlocal: 750 ave 760 max 735 min +Histogram: 1 0 0 0 1 0 0 0 0 2 +Nghost: 6230.5 ave 6253 max 6193 min +Histogram: 1 0 0 0 0 0 1 0 1 1 +Neighs: 276995 ave 280886 max 271360 min +Histogram: 1 0 0 0 1 0 0 0 1 1 + +Total # of neighbors = 1107981 +Ave neighs/atom = 369.327 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:03 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 new file mode 100644 index 0000000000..1187a755ee --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.1 @@ -0,0 +1,126 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.055 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 1 by 1 MPI processor grid + 3000 atoms + replicate CPU = 0.001 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +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} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +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 = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 539.2 | 539.2 | 539.2 Mbytes + Step Temp Press Density Volume + 0 300 10138.375 1 29915.273 + 10 295.97879 3575.2769 1 29915.273 + 20 292.76583 10309.128 1 29915.273 +Loop time of 10.8138 on 1 procs for 20 steps with 3000 atoms + +Performance: 0.080 ns/day, 300.383 hours/ns, 1.849 timesteps/s, 5.548 katom-step/s +99.9% CPU use with 1 MPI tasks x 1 OpenMP threads + +MPI task timing breakdown: +Section | min time | avg time | max time |%varavg| %total +--------------------------------------------------------------- +Pair | 4.7177 | 4.7177 | 4.7177 | 0.0 | 43.63 +Neigh | 0.17607 | 0.17607 | 0.17607 | 0.0 | 1.63 +Comm | 0.0017295 | 0.0017295 | 0.0017295 | 0.0 | 0.02 +Output | 8.5431e-05 | 8.5431e-05 | 8.5431e-05 | 0.0 | 0.00 +Modify | 5.9177 | 5.9177 | 5.9177 | 0.0 | 54.72 +Other | | 0.0004911 | | | 0.00 + +Nlocal: 3000 ave 3000 max 3000 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 11077 ave 11077 max 11077 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 971830 ave 971830 max 971830 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 971830 +Ave neighs/atom = 323.94333 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:12 diff --git a/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 new file mode 100644 index 0000000000..372156b6a2 --- /dev/null +++ b/examples/reaxff/water/log.29Aug24.reaxff.water-qtpie.g++.4 @@ -0,0 +1,126 @@ +LAMMPS (29 Aug 2024 - Development - patch_29Aug2024-76-g3f232caf9b) + using 1 OpenMP thread(s) per MPI task +# QTPIE Water + +boundary p p p +units real +atom_style charge + +read_data data.water +Reading data file ... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + reading atoms ... + 3000 atoms + read_data CPU = 0.053 seconds + +variable x index 1 +variable y index 1 +variable z index 1 + +replicate $x $y $z +replicate 1 $y $z +replicate 1 1 $z +replicate 1 1 1 +Replication is creating a 1x1x1 = 1 times larger system... + orthogonal box = (0 0 0) to (31.043046 31.043046 31.043046) + 1 by 2 by 2 MPI processor grid + 3000 atoms + replicate CPU = 0.002 seconds + +pair_style reaxff NULL safezone 3.0 mincap 150 +pair_coeff * * qeq_ff.water O H +WARNING: Changed valency_val to valency_boc for X (src/REAXFF/reaxff_ffield.cpp:294) +neighbor 0.5 bin +neigh_modify every 1 delay 0 check yes + +velocity all create 300.0 4928459 rot yes dist gaussian + +fix 1 all qtpie/reaxff 1 0.0 10.0 1.0e-6 reaxff gauss_exp.txt +fix 2 all nvt temp 300 300 50.0 + +timestep 0.5 +thermo 10 +thermo_style custom step temp press density vol + +run 20 + +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} +} + +- pair reaxff command: doi:10.1016/j.parco.2011.08.005 + +@Article{Aktulga12, + author = {H. M. Aktulga and J. C. Fogarty and S. A. Pandit and A. Y. Grama}, + title = {Parallel Reactive Molecular Dynamics: {N}umerical Methods and Algorithmic Techniques}, + journal = {Parallel Computing}, + year = 2012, + volume = 38, + number = {4--5}, + pages = {245--259} +} + +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 = 10.5 + ghost atom cutoff = 10.5 + binsize = 5.25, bins = 6 6 6 + 2 neighbor lists, perpetual/occasional/extra = 2 0 0 + (1) pair reaxff, perpetual + attributes: half, newton off, ghost + pair build: half/bin/ghost/newtoff + stencil: full/ghost/bin/3d + bin: standard + (2) fix qtpie/reaxff, perpetual, copy from (1) + attributes: half, newton off + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 260.5 | 262.2 | 263.6 Mbytes + Step Temp Press Density Volume + 0 300 10138.375 1 29915.273 + 10 295.97879 3575.2769 1 29915.273 + 20 292.76583 10309.128 1 29915.273 +Loop time of 3.13598 on 4 procs for 20 steps with 3000 atoms + +Performance: 0.276 ns/day, 87.111 hours/ns, 6.378 timesteps/s, 19.133 katom-step/s +99.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 | 1.6622 | 1.695 | 1.7252 | 2.2 | 54.05 +Neigh | 0.086543 | 0.087117 | 0.087848 | 0.2 | 2.78 +Comm | 0.0048192 | 0.035002 | 0.067754 | 15.4 | 1.12 +Output | 4.8033e-05 | 5.3375e-05 | 6.6893e-05 | 0.0 | 0.00 +Modify | 1.3176 | 1.3183 | 1.3189 | 0.0 | 42.04 +Other | | 0.0004753 | | | 0.02 + +Nlocal: 750 ave 760 max 735 min +Histogram: 1 0 0 0 1 0 0 0 0 2 +Nghost: 6229.5 ave 6253 max 6191 min +Histogram: 1 0 0 0 0 0 1 0 1 1 +Neighs: 277011 ave 280900 max 271380 min +Histogram: 1 0 0 0 1 0 0 0 1 1 + +Total # of neighbors = 1108044 +Ave neighs/atom = 369.348 +Neighbor list builds = 2 +Dangerous builds = 0 +Total wall time: 0:00:03 diff --git a/src/.gitignore b/src/.gitignore index b0c27849b2..f0554e3bfe 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -988,6 +988,8 @@ /fix_qeq_reaxff.h /fix_qmmm.cpp /fix_qmmm.h +/fix_qtpie_reaxff.cpp +/fix_qtpie_reaxff.h /fix_reaxff.cpp /fix_reaxff.h /fix_reaxff_bonds.cpp diff --git a/src/REAXFF/fix_qeq_reaxff.cpp b/src/REAXFF/fix_qeq_reaxff.cpp index 921f6e0261..067272f0f7 100644 --- a/src/REAXFF/fix_qeq_reaxff.cpp +++ b/src/REAXFF/fix_qeq_reaxff.cpp @@ -1158,7 +1158,7 @@ void FixQEqReaxFF::get_chi_field() for (int i = 0; i < nlocal; i++) { if (mask[i] & efgroupbit) { if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue; - chi_field[i] = -efield->efield[i][3]; + chi_field[i] = efield->efield[i][3]; } } } diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp new file mode 100644 index 0000000000..48c1109178 --- /dev/null +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -0,0 +1,1193 @@ +// 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 authors: + Efstratios M Kritikos, California Institute of Technology + (Implemented original version in LAMMMPS Aug 2019) + Navraj S Lalli, Imperial College London + (Reimplemented QTPIE as a new fix in LAMMPS Aug 2024 and extended functionality) +------------------------------------------------------------------------- */ + +#include "fix_qtpie_reaxff.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_efield.h" +#include "force.h" +#include "group.h" +#include "modify.h" +#include "neigh_list.h" +#include "neighbor.h" +#include "pair.h" +#include "region.h" +#include "respa.h" +#include "text_file_reader.h" +#include "update.h" + +#include "pair_reaxff.h" +#include "reaxff_api.h" + +#include +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +static constexpr double CONV_TO_EV = 14.4; +static constexpr double SMALL = 1.0e-14; +static constexpr double QSUMSMALL = 0.00001; +static constexpr double ANGSTROM_TO_BOHRRADIUS = 1.8897261259; + +/* ---------------------------------------------------------------------- */ + +FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), matvecs(0), pertype_option(nullptr), gauss_file(nullptr) +{ + // this fix returns a global scalar (the number of iterations) + scalar_flag = 1; + extscalar = 0; + + // this fix returns a per-atom vector (the effective electronegativity) + peratom_flag = 1; + size_peratom_cols = 0; + + imax = 200; + maxwarn = 1; + + if ((narg < 9) || (narg > 12)) error->all(FLERR,"Illegal fix {} command", style); + + nevery = utils::inumeric(FLERR,arg[3],false,lmp); + if (nevery <= 0) error->all(FLERR,"Illegal fix {} command", style); + + swa = utils::numeric(FLERR,arg[4],false,lmp); + swb = utils::numeric(FLERR,arg[5],false,lmp); + tolerance = utils::numeric(FLERR,arg[6],false,lmp); + pertype_option = utils::strdup(arg[7]); + gauss_file = utils::strdup(arg[8]); + + int iarg = 9; + while (iarg < narg) { + if (strcmp(arg[iarg],"nowarn") == 0) maxwarn = 0; + else if (strcmp(arg[iarg],"maxiter") == 0) { + if (iarg+1 > narg-1) + error->all(FLERR,"Illegal fix {} command", style); + imax = utils::numeric(FLERR,arg[iarg+1],false,lmp); + iarg++; + } else error->all(FLERR,"Illegal fix {} command", style); + iarg++; + } + shld = nullptr; + + nn = nt = n_cap = 0; + nmax = 0; + m_fill = m_cap = 0; + pack_flag = 0; + s = nullptr; + t = nullptr; + nprev = 4; + + Hdia_inv = nullptr; + b_s = nullptr; + chi_eff = nullptr; + b_t = nullptr; + b_prc = nullptr; + b_prm = nullptr; + + // CG + + p = nullptr; + q = nullptr; + r = nullptr; + d = nullptr; + + // H matrix + + H.firstnbr = nullptr; + H.numnbrs = nullptr; + H.jlist = nullptr; + H.val = nullptr; + + comm_forward = comm_reverse = 1; + + // perform initial allocation of atom-based arrays + // register with Atom class + + reaxff = dynamic_cast(force->pair_match("^reaxff",0)); + + s_hist = t_hist = nullptr; + atom->add_callback(Atom::GROW); +} + +/* ---------------------------------------------------------------------- */ + +FixQtpieReaxFF::~FixQtpieReaxFF() +{ + if (copymode) return; + + delete[] pertype_option; + delete[] gauss_file; + + // unregister callbacks to this fix from Atom class + + atom->delete_callback(id,Atom::GROW); + + memory->destroy(s_hist); + memory->destroy(t_hist); + + FixQtpieReaxFF::deallocate_storage(); + FixQtpieReaxFF::deallocate_matrix(); + + memory->destroy(shld); + + memory->destroy(gauss_exp); + if (!reaxflag) { + memory->destroy(chi); + memory->destroy(eta); + memory->destroy(gamma); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::post_constructor() +{ + grow_arrays(atom->nmax); + for (int i = 0; i < atom->nmax; i++) + for (int j = 0; j < nprev; ++j) + s_hist[i][j] = t_hist[i][j] = 0; + + pertype_parameters(pertype_option); +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::setmask() +{ + int mask = 0; + mask |= PRE_FORCE; + mask |= PRE_FORCE_RESPA; + mask |= MIN_PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pertype_parameters(char *arg) +{ + const int nlocal = atom->nlocal; + const int *mask = atom->mask; + const int *type = atom->type; + const int ntypes = atom->ntypes; + + // read gaussian orbital exponents + memory->create(gauss_exp,ntypes+1,"qtpie/reaxff:gauss_exp"); + if (comm->me == 0) { + gauss_exp[0] = 0.0; + try { + TextFileReader reader(gauss_file,"qtpie/reaxff gaussian exponents"); + reader.ignore_comments = true; + for (int i = 1; i <= ntypes; i++) { + const char *line = reader.next_line(); + if (!line) + throw TokenizerException("Fix qtpie/reaxff: Incorrect number of atom types in gauss file",""); + ValueTokenizer values(line); + + if (values.count() != 2) + throw TokenizerException("Fix qtpie/reaxff: Incorrect number of values per line " + "in gauss file",std::to_string(values.count())); + + int itype = values.next_int(); + if ((itype < 1) || (itype > ntypes)) + throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in gauss file", + std::to_string(itype)); + + double exp = values.next_double(); + if (exp < 0) + throw TokenizerException("Fix qtpie/reaxff: Invalid orbital exponent in gauss file", + std::to_string(exp)); + gauss_exp[itype] = exp; + } + } catch (std::exception &e) { + error->one(FLERR,e.what()); + } + } + + MPI_Bcast(gauss_exp,ntypes+1,MPI_DOUBLE,0,world); + + // define a cutoff distance (in atomic units) beyond which overlap integrals are neglected + // in calc_chi_eff() + const double exp_min = find_min_exp(gauss_exp,ntypes+1); + const int olap_cut = 10; // overlap integrals are neglected if less than pow(10,-olap_cut) + dist_cutoff = sqrt(2*olap_cut/exp_min*log(10.0)); + + // read chi, eta and gamma + + if (utils::strmatch(arg,"^reaxff")) { + reaxflag = 1; + Pair *pair = force->pair_match("^reaxff",0); + if (!pair) error->all(FLERR,"No reaxff pair style for fix qtpie/reaxff"); + + int tmp, tmp_all; + chi = (double *) pair->extract("chi",tmp); + eta = (double *) pair->extract("eta",tmp); + gamma = (double *) pair->extract("gamma",tmp); + if ((chi == nullptr) || (eta == nullptr) || (gamma == nullptr)) + error->all(FLERR, "Fix qtpie/reaxff could not extract qtpie parameters from pair reaxff"); + tmp = tmp_all = 0; + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + if ((chi[type[i]] == 0.0) && (eta[type[i]] == 0.0) && (gamma[type[i]] == 0.0)) + tmp = type[i]; + } + } + MPI_Allreduce(&tmp, &tmp_all, 1, MPI_INT, MPI_MAX, world); + if (tmp_all) + error->all(FLERR, "No qtpie parameters for atom type {} provided by pair reaxff", tmp_all); + return; + } else if (utils::strmatch(arg,"^reax/c")) { + error->all(FLERR, "Fix qtpie/reaxff keyword 'reax/c' is obsolete; please use 'reaxff'"); + } else if (platform::file_is_readable(arg)) { + ; // arg is readable file. will read below + } else { + error->all(FLERR, "Unknown fix qtpie/reaxff keyword {}", arg); + } + + reaxflag = 0; + + memory->create(chi,ntypes+1,"qtpie/reaxff:chi"); + memory->create(eta,ntypes+1,"qtpie/reaxff:eta"); + memory->create(gamma,ntypes+1,"qtpie/reaxff:gamma"); + + if (comm->me == 0) { + chi[0] = eta[0] = gamma[0] = 0.0; + try { + TextFileReader reader(arg,"qtpie/reaxff parameter"); + reader.ignore_comments = false; + for (int i = 1; i <= ntypes; i++) { + const char *line = reader.next_line(); + if (!line) + throw TokenizerException("Fix qtpie/reaxff: Invalid param file format",""); + ValueTokenizer values(line); + + if (values.count() != 4) + throw TokenizerException("Fix qtpie/reaxff: Incorrect format of param file",""); + + int itype = values.next_int(); + if ((itype < 1) || (itype > ntypes)) + throw TokenizerException("Fix qtpie/reaxff: Invalid atom type in param file", + std::to_string(itype)); + + chi[itype] = values.next_double(); + eta[itype] = values.next_double(); + gamma[itype] = values.next_double(); + } + } catch (std::exception &e) { + error->one(FLERR,e.what()); + } + } + + MPI_Bcast(chi,ntypes+1,MPI_DOUBLE,0,world); + MPI_Bcast(eta,ntypes+1,MPI_DOUBLE,0,world); + MPI_Bcast(gamma,ntypes+1,MPI_DOUBLE,0,world); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::allocate_storage() +{ + nmax = atom->nmax; + + memory->create(s,nmax,"qtpie:s"); + memory->create(t,nmax,"qtpie:t"); + + memory->create(Hdia_inv,nmax,"qtpie:Hdia_inv"); + memory->create(b_s,nmax,"qtpie:b_s"); + memory->create(chi_eff,nmax,"qtpie:chi_eff"); + vector_atom = chi_eff; + memory->create(b_t,nmax,"qtpie:b_t"); + memory->create(b_prc,nmax,"qtpie:b_prc"); + memory->create(b_prm,nmax,"qtpie:b_prm"); + + int size = nmax; + + memory->create(p,size,"qtpie:p"); + memory->create(q,size,"qtpie:q"); + memory->create(r,size,"qtpie:r"); + memory->create(d,size,"qtpie:d"); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::deallocate_storage() +{ + memory->destroy(s); + memory->destroy(t); + + memory->destroy(Hdia_inv); + memory->destroy(b_s); + memory->destroy(b_t); + memory->destroy(b_prc); + memory->destroy(b_prm); + memory->destroy(chi_eff); + + memory->destroy(p); + memory->destroy(q); + memory->destroy(r); + memory->destroy(d); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::reallocate_storage() +{ + deallocate_storage(); + allocate_storage(); + init_storage(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::allocate_matrix() +{ + int i,ii; + bigint m; + + int mincap; + double safezone; + + if (reaxflag) { + mincap = reaxff->api->system->mincap; + safezone = reaxff->api->system->safezone; + } else { + mincap = REAX_MIN_CAP; + safezone = REAX_SAFE_ZONE; + } + + n_cap = MAX((int)(atom->nlocal * safezone), mincap); + + // determine the total space for the H matrix + + m = 0; + for (ii = 0; ii < nn; ii++) { + i = ilist[ii]; + m += numneigh[i]; + } + bigint m_cap_big = (bigint)MAX(m * safezone, mincap * REAX_MIN_NBRS); + if (m_cap_big > MAXSMALLINT) + error->one(FLERR,"Too many neighbors in fix {}",style); + m_cap = m_cap_big; + + H.n = n_cap; + H.m = m_cap; + memory->create(H.firstnbr,n_cap,"qtpie:H.firstnbr"); + memory->create(H.numnbrs,n_cap,"qtpie:H.numnbrs"); + memory->create(H.jlist,m_cap,"qtpie:H.jlist"); + memory->create(H.val,m_cap,"qtpie:H.val"); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::deallocate_matrix() +{ + memory->destroy(H.firstnbr); + memory->destroy(H.numnbrs); + memory->destroy(H.jlist); + memory->destroy(H.val); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::reallocate_matrix() +{ + deallocate_matrix(); + allocate_matrix(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init() +{ + if (!atom->q_flag) + error->all(FLERR,"Fix {} requires atom attribute q", style); + + if (group->count(igroup) == 0) + error->all(FLERR,"Fix {} group has no atoms", style); + + // compute net charge and print warning if too large + double qsum_local = 0.0, qsum = 0.0; + for (int i = 0; i < atom->nlocal; i++) { + if (atom->mask[i] & groupbit) + qsum_local += atom->q[i]; + } + MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world); + + if ((comm->me == 0) && (fabs(qsum) > QSUMSMALL)) + error->warning(FLERR,"Fix {} group is not charge neutral, net charge = {:.8}", style, qsum); + + // get pointer to fix efield if present. there may be at most one instance of fix efield in use. + efield = nullptr; + auto fixes = modify->get_fix_by_style("^efield"); + if (fixes.size() == 1) efield = dynamic_cast(fixes.front()); + else if (fixes.size() > 1) + error->all(FLERR, "There may be only one fix efield instance used with fix {}", style); + + // ensure that fix efield is properly initialized before accessing its data and check some settings + if (efield) { + efield->init(); + if (strcmp(update->unit_style,"real") != 0) + error->all(FLERR,"Must use unit_style real with fix {} and external fields", style); + + if (efield->groupbit != 1){ // if efield is not applied to all atoms + error->all(FLERR,"Must use group id all for fix efield when using fix {}", style); + } + + if (efield->region){ // if efield is not applied to all atoms + error->all(FLERR,"Keyword region not supported for fix efield when using fix {}", style); + } + + if (efield->varflag == FixEfield::ATOM && efield->pstyle != FixEfield::ATOM) + error->all(FLERR,"Atom-style external electric field requires atom-style " + "potential variable when used with fix {}", style); + } + + // we need a half neighbor list w/ Newton off + // built whenever re-neighboring occurs + + neighbor->add_request(this, NeighConst::REQ_NEWTON_OFF); + + init_shielding(); + init_taper(); + + if (utils::strmatch(update->integrate_style,"^respa")) + nlevels_respa = (dynamic_cast(update->integrate))->nlevels; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::compute_scalar() +{ + return matvecs/2.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_shielding() +{ + int i,j; + int ntypes; + + ntypes = atom->ntypes; + if (shld == nullptr) + memory->create(shld,ntypes+1,ntypes+1,"qtpie:shielding"); + + for (i = 1; i <= ntypes; ++i) + for (j = 1; j <= ntypes; ++j) + shld[i][j] = pow(gamma[i] * gamma[j], -1.5); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_taper() +{ + double d7, swa2, swa3, swb2, swb3; + + if (fabs(swa) > 0.01 && comm->me == 0) + error->warning(FLERR,"Fix qtpie/reaxff has non-zero lower Taper radius cutoff"); + if (swb < 0) + error->all(FLERR, "Fix qtpie/reaxff has negative upper Taper radius cutoff"); + else if (swb < 5 && comm->me == 0) + error->warning(FLERR,"Fix qtpie/reaxff has very low Taper radius cutoff"); + + d7 = pow(swb - swa, 7); + swa2 = SQR(swa); + swa3 = CUBE(swa); + swb2 = SQR(swb); + swb3 = CUBE(swb); + + Tap[7] = 20.0 / d7; + Tap[6] = -70.0 * (swa + swb) / d7; + Tap[5] = 84.0 * (swa2 + 3.0*swa*swb + swb2) / d7; + Tap[4] = -35.0 * (swa3 + 9.0*swa2*swb + 9.0*swa*swb2 + swb3) / d7; + Tap[3] = 140.0 * (swa3*swb + 3.0*swa2*swb2 + swa*swb3) / d7; + Tap[2] =-210.0 * (swa3*swb2 + swa2*swb3) / d7; + Tap[1] = 140.0 * swa3 * swb3 / d7; + Tap[0] = (-35.0*swa3*swb2*swb2 + 21.0*swa2*swb3*swb2 - + 7.0*swa*swb3*swb3 + swb3*swb3*swb) / d7; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::setup_pre_force(int vflag) +{ + if (reaxff) { + nn = reaxff->list->inum; + nt = reaxff->list->inum + reaxff->list->gnum; + ilist = reaxff->list->ilist; + numneigh = reaxff->list->numneigh; + firstneigh = reaxff->list->firstneigh; + } else { + nn = list->inum; + nt = list->inum + list->gnum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + } + + deallocate_storage(); + allocate_storage(); + + init_storage(); + + deallocate_matrix(); + allocate_matrix(); + + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::setup_pre_force_respa(int vflag, int ilevel) +{ + if (ilevel < nlevels_respa-1) return; + setup_pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::min_setup_pre_force(int vflag) +{ + setup_pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_storage() +{ + calc_chi_eff(); + + for (int ii = 0; ii < nn; ii++) { + int i = ilist[ii]; + if (atom->mask[i] & groupbit) { + Hdia_inv[i] = 1. / eta[atom->type[i]]; + b_s[i] = -chi_eff[i]; + b_t[i] = -1.0; + b_prc[i] = 0; + b_prm[i] = 0; + s[i] = t[i] = 0; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pre_force(int /*vflag*/) +{ + if (update->ntimestep % nevery) return; + + int n = atom->nlocal; + + if (reaxff) { + nn = reaxff->list->inum; + nt = reaxff->list->inum + reaxff->list->gnum; + ilist = reaxff->list->ilist; + numneigh = reaxff->list->numneigh; + firstneigh = reaxff->list->firstneigh; + } else { + nn = list->inum; + nt = list->inum + list->gnum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + } + + // grow arrays if necessary + // need to be atom->nmax in length + + if (atom->nmax > nmax) reallocate_storage(); + if (n > n_cap*DANGER_ZONE || m_fill > m_cap*DANGER_ZONE) + reallocate_matrix(); + + calc_chi_eff(); + + init_matvec(); + + matvecs_s = CG(b_s, s); // CG on s - parallel + matvecs_t = CG(b_t, t); // CG on t - parallel + matvecs = matvecs_s + matvecs_t; + + calculate_Q(); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::pre_force_respa(int vflag, int ilevel, int /*iloop*/) +{ + if (ilevel == nlevels_respa-1) pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::min_pre_force(int vflag) +{ + pre_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::init_matvec() +{ + /* fill-in H matrix */ + compute_H(); + + int ii, i; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + + /* init pre-conditioner for H and init solution vectors */ + Hdia_inv[i] = 1. / eta[atom->type[i]]; + b_s[i] = -chi_eff[i]; + b_t[i] = -1.0; + + /* quadratic extrapolation for s & t from previous solutions */ + t[i] = t_hist[i][2] + 3 * (t_hist[i][0] - t_hist[i][1]); + + /* cubic extrapolation for s & t from previous solutions */ + 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 FixQtpieReaxFF::compute_H() +{ + int jnum; + int i, j, ii, jj, flag; + double dx, dy, dz, r_sqr; + constexpr double EPSILON = 0.0001; + + int *type = atom->type; + tagint *tag = atom->tag; + double **x = atom->x; + int *mask = atom->mask; + + // fill in the H matrix + m_fill = 0; + r_sqr = 0; + for (ii = 0; ii < nn; 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 = SQR(dx) + SQR(dy) + SQR(dz); + + flag = 0; + if (r_sqr <= SQR(swb)) { + if (j < atom->nlocal) flag = 1; + else if (tag[i] < tag[j]) flag = 1; + else if (tag[i] == tag[j]) { + if (dz > EPSILON) flag = 1; + else if (fabs(dz) < EPSILON) { + if (dy > EPSILON) flag = 1; + else if (fabs(dy) < EPSILON && dx > EPSILON) + flag = 1; + } + } + } + + if (flag) { + H.jlist[m_fill] = j; + H.val[m_fill] = calculate_H(sqrt(r_sqr), shld[type[i]][type[j]]); + m_fill++; + } + } + H.numnbrs[i] = m_fill - H.firstnbr[i]; + } + } + + if (m_fill >= H.m) + error->all(FLERR,"Fix qtpie/reaxff H matrix size has been exceeded: m_fill={} H.m={}\n", + m_fill, H.m); +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::calculate_H(double r, double gamma) +{ + double Taper, denom; + + Taper = Tap[7] * r + Tap[6]; + Taper = Taper * r + Tap[5]; + Taper = Taper * r + Tap[4]; + Taper = Taper * r + Tap[3]; + Taper = Taper * r + Tap[2]; + Taper = Taper * r + Tap[1]; + Taper = Taper * r + Tap[0]; + + denom = r * r * r + gamma; + denom = pow(denom,1.0/3.0); + + return Taper * CONV_TO_EV / denom; +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::CG(double *b, double *x) +{ + int i, j; + double tmp, alpha, beta, b_norm; + double sig_old, sig_new; + + int jj; + + pack_flag = 1; + sparse_matvec(&H, x, q); + comm->reverse_comm(this); //Coll_Vector(q); + + vector_sum(r , 1., b, -1., q, nn); + + for (jj = 0; jj < nn; ++jj) { + j = ilist[jj]; + if (atom->mask[j] & groupbit) + d[j] = r[j] * Hdia_inv[j]; //pre-condition + } + + b_norm = parallel_norm(b, nn); + sig_new = parallel_dot(r, d, nn); + + for (i = 1; i < imax && sqrt(sig_new) / b_norm > tolerance; ++i) { + comm->forward_comm(this); //Dist_vector(d); + sparse_matvec(&H, d, q); + comm->reverse_comm(this); //Coll_vector(q); + + tmp = parallel_dot(d, q, nn); + alpha = sig_new / tmp; + + vector_add(x, alpha, d, nn); + vector_add(r, -alpha, q, nn); + + // pre-conditioning + for (jj = 0; jj < nn; ++jj) { + j = ilist[jj]; + if (atom->mask[j] & groupbit) + p[j] = r[j] * Hdia_inv[j]; + } + + sig_old = sig_new; + sig_new = parallel_dot(r, p, nn); + + beta = sig_new / sig_old; + vector_sum(d, 1., p, beta, d, nn); + } + + if ((i >= imax) && maxwarn && (comm->me == 0)) + error->warning(FLERR, "Fix qtpie/reaxff CG convergence failed after {} iterations at step {}", + i,update->ntimestep); + return i; +} + + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::sparse_matvec(sparse_matrix *A, double *x, double *b) +{ + int i, j, itr_j; + int ii; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + b[i] = eta[atom->type[i]] * x[i]; + } + + int nall = atom->nlocal + atom->nghost; + for (i = atom->nlocal; i < nall; ++i) + b[i] = 0; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + for (itr_j=A->firstnbr[i]; itr_jfirstnbr[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]; + } + } + } + +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::calculate_Q() +{ + int i, k; + double u, s_sum, t_sum; + double *q = atom->q; + + int ii; + + s_sum = parallel_vector_acc(s, nn); + t_sum = parallel_vector_acc(t, nn); + u = s_sum / t_sum; + + for (ii = 0; ii < nn; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) { + q[i] = s[i] - u * t[i]; + + /* backup s & t */ + for (k = nprev-1; 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]; + } + } + + pack_flag = 4; + comm->forward_comm(this); //Dist_vector(atom->q); +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int m; + + if (pack_flag == 1) + for (m = 0; m < n; m++) buf[m] = d[list[m]]; + else if (pack_flag == 2) + for (m = 0; m < n; m++) buf[m] = s[list[m]]; + else if (pack_flag == 3) + for (m = 0; m < n; m++) buf[m] = t[list[m]]; + else if (pack_flag == 4) + for (m = 0; m < n; m++) buf[m] = atom->q[list[m]]; + return n; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m; + + if (pack_flag == 1) + for (m = 0, i = first; m < n; m++, i++) d[i] = buf[m]; + else if (pack_flag == 2) + for (m = 0, i = first; m < n; m++, i++) s[i] = buf[m]; + else if (pack_flag == 3) + for (m = 0, i = first; m < n; m++, i++) t[i] = buf[m]; + else if (pack_flag == 4) + for (m = 0, i = first; m < n; m++, i++) atom->q[i] = buf[m]; +} + +/* ---------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m; + for (m = 0, i = first; m < n; m++, i++) buf[m] = q[i]; + return n; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::unpack_reverse_comm(int n, int *list, double *buf) +{ + for (int m = 0; m < n; m++) q[list[m]] += buf[m]; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double FixQtpieReaxFF::memory_usage() +{ + double bytes; + + bytes = (double)atom->nmax*nprev*2 * sizeof(double); // s_hist & t_hist + bytes += (double)atom->nmax*11 * sizeof(double); // storage + bytes += (double)n_cap*2 * sizeof(int); // matrix... + bytes += (double)m_cap * sizeof(int); + bytes += (double)m_cap * sizeof(double); + + return bytes; +} + +/* ---------------------------------------------------------------------- + allocate fictitious charge arrays +------------------------------------------------------------------------- */ + +void FixQtpieReaxFF::grow_arrays(int nmax) +{ + memory->grow(s_hist,nmax,nprev,"qtpie:s_hist"); + memory->grow(t_hist,nmax,nprev,"qtpie:t_hist"); +} + +/* ---------------------------------------------------------------------- + copy values within fictitious charge arrays +------------------------------------------------------------------------- */ + +void FixQtpieReaxFF::copy_arrays(int i, int j, int /*delflag*/) +{ + for (int m = 0; m < nprev; m++) { + s_hist[j][m] = s_hist[i][m]; + t_hist[j][m] = t_hist[i][m]; + } +} + +/* ---------------------------------------------------------------------- + pack values in local atom-based array for exchange with another proc +------------------------------------------------------------------------- */ + +int FixQtpieReaxFF::pack_exchange(int i, double *buf) +{ + for (int m = 0; m < nprev; m++) buf[m] = s_hist[i][m]; + for (int m = 0; m < nprev; m++) buf[nprev+m] = t_hist[i][m]; + return nprev*2; +} + +/* ---------------------------------------------------------------------- + unpack values in local atom-based array from exchange with another proc +------------------------------------------------------------------------- */ + +int FixQtpieReaxFF::unpack_exchange(int nlocal, double *buf) +{ + for (int m = 0; m < nprev; m++) s_hist[nlocal][m] = buf[m]; + for (int m = 0; m < nprev; m++) t_hist[nlocal][m] = buf[nprev+m]; + return nprev*2; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_norm(double *v, int n) +{ + int i; + double my_sum, norm_sqr; + + int ii; + + my_sum = 0.0; + norm_sqr = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_sum += SQR(v[i]); + } + + MPI_Allreduce(&my_sum, &norm_sqr, 1, MPI_DOUBLE, MPI_SUM, world); + + return sqrt(norm_sqr); +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_dot(double *v1, double *v2, int n) +{ + int i; + double my_dot, res; + + int ii; + + my_dot = 0.0; + res = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_dot += v1[i] * v2[i]; + } + + MPI_Allreduce(&my_dot, &res, 1, MPI_DOUBLE, MPI_SUM, world); + + return res; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::parallel_vector_acc(double *v, int n) +{ + int i; + double my_acc, res; + + int ii; + + my_acc = 0.0; + res = 0.0; + for (ii = 0; ii < n; ++ii) { + i = ilist[ii]; + if (atom->mask[i] & groupbit) + my_acc += v[i]; + } + + MPI_Allreduce(&my_acc, &res, 1, MPI_DOUBLE, MPI_SUM, world); + + return res; +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::vector_sum(double* dest, double c, double* v, + double d, double* y, int k) +{ + int kk; + + for (--k; k>=0; --k) { + kk = ilist[k]; + if (atom->mask[kk] & groupbit) + dest[kk] = c * v[kk] + d * y[kk]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::vector_add(double* dest, double c, double* v, int k) +{ + int kk; + + for (--k; k>=0; --k) { + kk = ilist[k]; + if (atom->mask[kk] & groupbit) + dest[kk] += c * v[kk]; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixQtpieReaxFF::calc_chi_eff() +{ + memset(&chi_eff[0],0,atom->nmax*sizeof(double)); + + const auto x = (const double * const *)atom->x; + const int ntypes = atom->ntypes; + const int *type = atom->type; + + double dist,overlap,sum_n,sum_d,expa,expb,chia,chib,phia,phib,p,m; + int i,j; + + // check ghost atoms are stored up to the distance cutoff for overlap integrals + const double comm_cutoff = MAX(neighbor->cutneighmax,comm->cutghostuser); + if(comm_cutoff < dist_cutoff/ANGSTROM_TO_BOHRRADIUS) { + error->all(FLERR,"comm cutoff = {} Angstrom is smaller than distance cutoff = {} Angstrom " + "for overlap integrals in {}. Increase comm cutoff with comm_modify", + comm_cutoff, dist_cutoff/ANGSTROM_TO_BOHRRADIUS, style); + } + + // efield energy is in real units of kcal/mol, factor needed for conversion to eV + const double qe2f = force->qe2f; + const double factor = 1.0/qe2f; + + if (efield) { + if (efield->varflag != FixEfield::CONSTANT) + efield->update_efield_variables(); + } + + // compute chi_eff for each local atom + for (i = 0; i < nn; i++) { + expa = gauss_exp[type[i]]; + chia = chi[type[i]]; + if (efield) { + if (efield->varflag != FixEfield::ATOM) { + phia = -factor*(x[i][0]*efield->ex + x[i][1]*efield->ey + x[i][2]*efield->ez); + } else { // atom-style potential from FixEfield + phia = efield->efield[i][3]; + } + } + + sum_n = 0.0; + sum_d = 0.0; + + for (j = 0; j < nt; j++) { + dist = distance(x[i],x[j])*ANGSTROM_TO_BOHRRADIUS; // in atomic units + + if (dist < dist_cutoff) { + expb = gauss_exp[type[j]]; + chib = chi[type[j]]; + + // overlap integral of two normalised 1s Gaussian type orbitals + p = expa + expb; + m = expa * expb / p; + overlap = pow((4.0*m/p),0.75) * exp(-m*dist*dist); + + if (efield) { + if (efield->varflag != FixEfield::ATOM) { + phib = -factor*(x[j][0]*efield->ex + x[j][1]*efield->ey + x[j][2]*efield->ez); + } else { // atom-style potential from FixEfield + phib = efield->efield[j][3]; + } + sum_n += (chia - chib + phia - phib) * overlap; + } else { + sum_n += (chia - chib) * overlap; + } + sum_d += overlap; + } + } + + chi_eff[i] = sum_n / sum_d; + } +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::find_min_exp(const double *array, const int array_length) +{ + // index of first gaussian orbital exponent is 1 + double exp_min = array[1]; + for (int i = 2; i < array_length; i++) + { + if (array[i] < exp_min) + exp_min = array[i]; + } + return exp_min; +} + +/* ---------------------------------------------------------------------- */ + +double FixQtpieReaxFF::distance(const double *posa, const double *posb) +{ + double dx, dy, dz; + dx = posb[0] - posa[0]; + dy = posb[1] - posa[1]; + dz = posb[2] - posa[2]; + return sqrt(dx*dx + dy*dy + dz*dz); +} diff --git a/src/REAXFF/fix_qtpie_reaxff.h b/src/REAXFF/fix_qtpie_reaxff.h new file mode 100644 index 0000000000..2d82ad197c --- /dev/null +++ b/src/REAXFF/fix_qtpie_reaxff.h @@ -0,0 +1,141 @@ +/* -*- 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(qtpie/reaxff,FixQtpieReaxFF); +// clang-format on +#else + +#ifndef LMP_FIX_QTPIE_REAXFF_H +#define LMP_FIX_QTPIE_REAXFF_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixQtpieReaxFF : public Fix { + public: + FixQtpieReaxFF(class LAMMPS *, int, char **); + ~FixQtpieReaxFF() override; + int setmask() override; + void post_constructor() override; + void init() override; + void init_list(int, class NeighList *) override; + void init_storage(); + void setup_pre_force(int) override; + void pre_force(int) override; + + void setup_pre_force_respa(int, int) override; + void pre_force_respa(int, int, int) override; + + void min_setup_pre_force(int); + void min_pre_force(int) override; + + double compute_scalar() override; + + protected: + int nevery, reaxflag; + int matvecs; + int nn, nt, m_fill; + int n_cap, nmax, m_cap; + int pack_flag; + int nlevels_respa; + class NeighList *list; + class PairReaxFF *reaxff; + class FixEfield *efield; + int *ilist, *jlist, *numneigh, **firstneigh; + + double swa, swb; // lower/upper Taper cutoff radius + double Tap[8]; // Taper function + double tolerance; // tolerance for the norm of the rel residual in CG + + double *chi, *eta, *gamma; // qeq parameters + double **shld; + + // fictitious charges + + double *s, *t; + double **s_hist, **t_hist; + int nprev; + + typedef struct { + int n, m; + int *firstnbr; + int *numnbrs; + int *jlist; + double *val; + } sparse_matrix; + + sparse_matrix H; + double *Hdia_inv; + double *b_s, *b_t; + double *b_prc, *b_prm; + double *chi_eff; // array of effective electronegativities + + //CG storage + double *p, *q, *r, *d; + int imax, maxwarn; + + char *pertype_option; // argument to determine how per-type info is obtained + char *gauss_file; // input file for gaussian orbital exponents + double *gauss_exp; // array of gaussian orbital exponents for each atom type + double dist_cutoff; // separation distance beyond which to neglect overlap integrals + + void pertype_parameters(char *); + void init_shielding(); + void init_taper(); + void allocate_storage(); + void deallocate_storage(); + void reallocate_storage(); + void allocate_matrix(); + void deallocate_matrix(); + void reallocate_matrix(); + + void init_matvec(); + void init_H(); + void compute_H(); + double calculate_H(double, double); + void calculate_Q(); + + int CG(double *, double *); + void sparse_matvec(sparse_matrix *, double *, double *); + + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + double memory_usage() override; + void grow_arrays(int) override; + void copy_arrays(int, int, int) override; + int pack_exchange(int, double *) override; + int unpack_exchange(int, double *) override; + + double parallel_norm(double *, int); + double parallel_dot(double *, double *, int); + double parallel_vector_acc(double *, int); + + void vector_sum(double *, double, double *, double, double *, int); + void vector_add(double *, double, double *, int); + + void calc_chi_eff(); + double find_min_exp(const double*, const int); + double distance(const double*, const double*); + + int matvecs_s, matvecs_t; // Iteration count for each system +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/REAXFF/pair_reaxff.cpp b/src/REAXFF/pair_reaxff.cpp index 08e90933b2..ca60e5248f 100644 --- a/src/REAXFF/pair_reaxff.cpp +++ b/src/REAXFF/pair_reaxff.cpp @@ -340,11 +340,13 @@ void PairReaxFF::init_style() auto acks2_fixes = modify->get_fix_by_style("^acks2/reax"); int have_qeq = modify->get_fix_by_style("^qeq/reax").size() - + modify->get_fix_by_style("^qeq/shielded").size() + acks2_fixes.size(); + + modify->get_fix_by_style("^qeq/shielded").size() + acks2_fixes.size() + + modify->get_fix_by_style("^qtpie/reax").size(); if (qeqflag && (have_qeq != 1)) error->all(FLERR,"Pair style reaxff requires use of exactly one of the " - "fix qeq/reaxff or fix qeq/shielded or fix acks2/reaxff commands"); + "fix qeq/reaxff or fix qeq/shielded or fix acks2/reaxff or " + "fix qtpie/reaxff commands"); api->system->acks2_flag = acks2_fixes.size(); if (api->system->acks2_flag) diff --git a/src/fix_efield.h b/src/fix_efield.h index 72fd204898..108395cc2c 100644 --- a/src/fix_efield.h +++ b/src/fix_efield.h @@ -26,6 +26,7 @@ namespace LAMMPS_NS { class FixEfield : public Fix { friend class FixQEqReaxFF; + friend class FixQtpieReaxFF; public: FixEfield(class LAMMPS *, int, char **);