doc pages and sync qm vs qmmm fixes

This commit is contained in:
Steve Plimpton
2023-02-22 17:21:40 -07:00
parent 17e1ed4869
commit 8b16301e5f
12 changed files with 5308 additions and 667 deletions

View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`fix <fix>` command
* mdi/qm = style name of this fix command
* zero or more keyword/value pairs may be appended
* keyword = *virial* or *add* or *every* or *connect* or *elements*
* keyword = *virial* or *add* or *every* or *connect* or *elements* or *mc*
.. parsed-literal::
@ -29,7 +29,9 @@ Syntax
yes = perform a one-time connection to the MDI engine code
no = do not perform the connection operation
*elements* args = N_1 N_2 ... N_ntypes
N_1,N_2,...N_ntypes = atomic number for each of ntypes LAMMPS atom types
N_1,N_2,...N_ntypes = chemical symbol for each of ntypes LAMMPS atom types
*mc* args = mcfixID
mcfixID = ID of a Monte Carlo fix designed to work with this fix
Examples
""""""""
@ -38,7 +40,7 @@ Examples
fix 1 all mdi/qm
fix 1 all mdi/qm virial yes
fix 1 all mdi/qm add no every 100 elements 13 29
fix 1 all mdi/qm add no every 100 elements C C H O
Description
"""""""""""
@ -57,12 +59,25 @@ The server code must support use of the `MDI Library
<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ as
explained below.
Typically, to use this fix, the input script should not define any
other classical force field components, e.g. a pair style, bond style,
etc.
These are example use cases for this fix, discussed further below:
* perform an ab initio MD (AIMD) simulation with quantum forces
* perform an energy minimization with quantum forces
* perform a nudged elastic band (NEB) calculation with quantum forces
* perform a QM calculation for a series of independent systems which LAMMPS reads or generates
* perform a QM calculation for a series of independent systems which LAMMPS reads or generates once
* run a classical MD simulation and calculate QM energy/forces once every N steps on the current configuration
More generally any command which calulates per-atom forces can instead
use quantum forces be defining this fix. Examples are the Monte Carlo
commands :doc:`fix gcmc <fix_gcmc>` and :doc:`fix atom/swap
<fix_atom_swap>`, as well as the :doc:`compute born <compute_born>`
command. The only requirement is that internally the commmand invokes
the post_force() method of fixes such as this one, which will trigger
the quantum calculation.
The code coupling performed by this command is done via the `MDI
Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
@ -80,19 +95,23 @@ to launch two codes so that they communicate via the MDI library using
either MPI or sockets. Any QM code that supports MDI could be used in
place of LAMMPS acting as a QM surrogate. See the :doc:`Howto mdi
<Howto_mdi>` page for a current list (March 2022) of such QM codes.
The examples/QUANTUM directory has examples for coupling LAMMPS to 3
QM codes either via this fix or the :doc:`fix mdi/qmmm <fix_mdi_qmmm>`
command.
Note that an engine code can support MDI in either or both of two
modes. It can be used as a stand-alone code, launched at the same
time as LAMMPS. Or it can be used as a plugin library, which LAMMPS
loads. See the :doc:`mdi plugin <mdi>` command for how to trigger
LAMMPS to load a plugin library. The examples/mdi/README file
explains how to launch the two codes in either mode.
LAMMPS to load a plugin library. The examples/mdi/README file and
examples/QUANTUM/QM-code/README files explain how to launch the two
codes in either mode.
----------
The *virial* keyword setting of yes or no determines whether
LAMMPS will request the QM code to also compute and return
a 6-element symmetric virial tensor for the system.
The *virial* keyword setting of yes or no determines whether LAMMPS
will request the QM code to also compute and return a symmetric virial
tensor for the system.
The *add* keyword setting of *yes* or *no* determines whether the
energy and forces and virial returned by the QM code will be added to
@ -109,25 +128,27 @@ commands. See details below.
The *every* keyword determines how often the QM code will be invoked
during a dynamics run with the current LAMMPS simulation box and
configuration of atoms. The QM code will be called once every
*Nevery* timesteps.
*Nevery* timesteps. By default *Nevery* = 1.
The *connect* keyword determines whether this fix performs a one-time
connection to the QM code. The default is *yes*. The only time a
*no* is needed is if this command is used multiple times in an input
script. E.g. if it used inside a loop which also uses the :doc:`clear
<clear>` command to destroy the system (including any defined fixes).
See the examples/mdi/in.series.driver script as an example of this,
where LAMMPS is using the QM code to compute energy and forces for a
series of system configurations. In this use case *connect no*
is used along with the :doc:`mdi connect and exit <mdi>` command
to one-time initiate/terminate the connection outside the loop.
script and the MDI coupling is between two stand-alone codes (not
plugin mode). E.g. if it used inside a loop which also uses the
:doc:`clear <clear>` command to destroy the system (including this
fix). See the examples/mdi/in.series.driver script as an example of
this, where LAMMPS is using the QM code to compute energy and forces
for a series of system configurations. In this use case *connect no*
is used along with the :doc:`mdi connect and exit <mdi>` command to
one-time initiate/terminate the connection outside the loop.
The *elements* keyword allows specification of what element each
LAMMPS atom type corresponds to. This is specified by the atomic
number of the element, e.g. 13 for Al. An atomic number must be
specified for each of the ntypes LAMMPS atom types. Ntypes is
typically specified via the create_box command or in the data file
read by the read_data command.
LAMMPS atom type corresponds to. This is specified by the chemical
symbol of the element, e.g. C or Al or Si. A symbol must be specified
for each of the ntypes LAMMPS atom types. Multiple LAMMPS types can
represent the same element. Ntypes is typically specified via the
:doc:`create_box <create_box>` command or in the data file read by the
:doc:`read_data <read_data>` command.
If this keyword is specified, then this fix will send the MDI
">ELEMENTS" command to the engine, to ensure the two codes are
@ -136,6 +157,14 @@ not specified, then this fix will send the MDI >TYPES command to the
engine. This is fine if both the LAMMPS driver and the MDI engine are
initialized so that the atom type values are consistent in both codes.
The *mc* keyword enables this fix to be used with a Monte Carlo (MC)
fix to calculate before/after quantum energies as part of the MC
accept/reject criterion. The :doc:`fix gcmc <fix_gcmc>` and :doc:`fix
atom/swap <fix_atom_swap>` commands can be used in this manner.
Specify the ID of the MC fix following the *mc* keyword. This allows
the two fixes to coordinate when MC events are being calculated versus
MD timesteps between the MC events.
----------
The following 3 example use cases are illustrated in the examples/mdi
@ -265,7 +294,8 @@ unit conversions between LAMMPS and MDI units. The other code will
also perform similar unit conversions into its preferred units.
LAMMPS can also be used as an MDI driver in other unit choices it
supports, e.g. *lj*, but then no unit conversion is performed.
supports, e.g. *lj*, but then no unit conversion to MDI units is
performed.
Related commands
""""""""""""""""

View File

@ -1,6 +1,6 @@
.. index:: fix mdi/qm
.. index:: fix mdi/qmmm
fix mdi/qm command
fix mdi/qmmm command
======================
Syntax
@ -11,7 +11,7 @@ Syntax
fix ID group-ID mdi/qmmm mode keyword value(s) keyword value(s) ...
* ID, group-ID are documented in :doc:`fix <fix>` command
* mdi/qm = style name of this fix command
* mdi/qmmm = style name of this fix command
* mode = *direct* or *potential*
* zero or more keyword/value pairs may be appended
* keyword = *virial* or *add* or *every* or *connect* or *elements*
@ -25,7 +25,7 @@ Syntax
yes = perform a one-time connection to the MDI engine code
no = do not perform the connection operation
*elements* args = N_1 N_2 ... N_ntypes
N_1,N_2,...N_ntypes = atomic number for each of ntypes LAMMPS atom types
N_1,N_2,...N_ntypes = chemical symbol for each of ntypes LAMMPS atom types
Examples
""""""""
@ -40,14 +40,29 @@ Description
"""""""""""
This command enables LAMMPS to act as a client with another server
code to perform a couple QMMM (quantum-mechanics molecular-mechanics)
simulation. LAMMPS will perform classical MD (molecular-mechnanics)
for the (typically larger) MM portion of the system. A quantum
mechanics code will calculate quantum effects for the QM portion of
the system. The QM server code must support use of the `MDI Library
code to perform a coupled QMMM (quantum-mechanics/molecular-mechanics)
simulation. LAMMPS will perform classical MD (molecular-mechnanics
or MM) for the (typically larger) MM portion of the system. A quantum
mechanics code will calculate quantum energy and forces for the QM
portion of the system. The two codes work together to calculate the
energy and forces due to the cross interactions between QM and MM
atoms. The QM server code must support use of the `MDI Library
<https://molssi-mdi.github.io/MDI_Library/html/index.html>`_ as
explained below.
The partitioning of the system between QM and MM atoms is determined
by the specified group. Atoms in that group are QM atoms; the remaining
atoms are MM atoms. The input script should
need to remove bonds/etc between QM atoms
currently no bonds between QM/MM atoms are allowed
need to be able to compute Coulomb only portion of the force
field in a pair style, use hybrid
must be a charged system
explain the 2 algs: DIRECT and POTENTIAL
The code coupling performed by this command is done via the `MDI
Library <https://molssi-mdi.github.io/MDI_Library/html/index.html>`_.
LAMMPS runs as an MDI driver (client), and sends MDI commands to an
@ -56,11 +71,6 @@ support for MDI. See the :doc:`Howto mdi <Howto_mdi>` page for more
information about how LAMMPS can operate as either an MDI driver or
engine.
Q: provide an example where LAMMPS is also the QM code ?
similar to fix mdi/qm ?
The examples/mdi directory contains input scripts using this fix in
the various use cases discussed below. In each case, two instances of
LAMMPS are used, once as an MDI driver, once as an MDI engine
@ -79,6 +89,8 @@ explains how to launch the two codes in either mode.
----------
----------
The *virial* keyword setting of yes or no determines whether
LAMMPS will request the QM code to also compute and return
a 6-element symmetric virial tensor for the system.
@ -87,19 +99,20 @@ The *connect* keyword determines whether this fix performs a one-time
connection to the QM code. The default is *yes*. The only time a
*no* is needed is if this command is used multiple times in an input
script. E.g. if it used inside a loop which also uses the :doc:`clear
<clear>` command to destroy the system (including any defined fixes).
See the examples/mdi/in.series.driver script as an example of this,
where LAMMPS is using the QM code to compute energy and forces for a
series of system configurations. In this use case *connect no*
is used along with the :doc:`mdi connect and exit <mdi>` command
to one-time initiate/terminate the connection outside the loop.
<clear>` command to destroy the system (including this fix). As
example would be a script which loop over a series of independent QMMM
simulations, e.g. each with their own data file. In this use case
*connect no* could be used along with the :doc:`mdi connect and exit
<mdi>` command to one-time initiate/terminate the connection outside
the loop.
The *elements* keyword allows specification of what element each
LAMMPS atom type corresponds to. This is specified by the atomic
number of the element, e.g. 13 for Al. An atomic number must be
specified for each of the ntypes LAMMPS atom types. Ntypes is
typically specified via the create_box command or in the data file
read by the read_data command.
LAMMPS atom type corresponds to. This is specified by the chemical
symbol of the element, e.g. C or Al or Si. A symbol must be specified
for each of the ntypes LAMMPS atom types. Multiple LAMMPS types can
represent the same element. Ntypes is typically specified via the
:doc:`create_box <create_box>` command or in the data file read by the
:doc:`read_data <read_data>` command.
If this keyword is specified, then this fix will send the MDI
">ELEMENTS" command to the engine, to insure the two codes are
@ -236,9 +249,6 @@ used to specify *real* or *metal* units. This will ensure the correct
unit conversions between LAMMPS and MDI units. The other code will
also perform similar unit conversions into its preferred units.
LAMMPS can also be used as an MDI driver in other unit choices it
supports, e.g. *lj*, but then no unit conversion is performed.
Related commands
""""""""""""""""

View File

@ -56,4 +56,4 @@ thermo_style custom step cpu temp ke evdwl ecoul epair emol elong &
f_2 pe etotal press
thermo 1
run 2
run 10

View File

@ -60,4 +60,4 @@ variable p equal extract_setting(world_size)
mdi plugin nwchem_mdi mdi "-role ENGINE -name NWChem -method LINK" &
extra "template.water.nw water.nw log.water.pwdft.qmmm.plugin.$p" &
command "run 1"
command "run 10"

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -94,7 +94,7 @@ fix_modify 2 energy yes
thermo_style custom step cpu temp ke evdwl ecoul epair emol elong f_2 pe etotal press
thermo 1
run 2
run 10
Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
@ -115,24 +115,32 @@ Neighbor list info ...
bin: standard
Per MPI rank memory allocation (min/avg/max) = 7.445 | 7.445 | 7.445 Mbytes
Step CPU Temp KinEng E_vdwl E_coul E_pair E_mol E_long f_2 PotEng TotEng Press
0 0 0 0 0.95937448 -170.68735 -169.72797 3.3962859 0 -10627.878 -10794.209 -10794.209 -374.81058
1 1.3827511 1.1713705 0.017458165 0.95934217 -170.68987 -169.73053 3.3804649 0 -10627.878 -10794.229 -10794.211 -376.91596
2 4.3658207 4.6633575 0.069502914 0.95924532 -170.70698 -169.74773 3.3333254 0 -10627.871 -10794.286 -10794.216 -383.31463
Loop time of 4.36586 on 1 procs for 2 steps with 6 atoms
0 0 0 0 0.95937448 -170.68735 -169.72798 3.3962859 0 -10627.878 -10794.209 -10794.209 -374.81058
1 1.3998539 1.1713705 0.017458165 0.95934217 -170.68987 -169.73053 3.3804649 0 -10627.878 -10794.229 -10794.211 -376.91596
2 4.2727994 4.6633575 0.069502914 0.95924532 -170.70698 -169.74773 3.3333254 0 -10627.871 -10794.286 -10794.216 -383.31463
3 8.6537003 10.410148 0.15515337 0.95908482 -170.73725 -169.77816 3.2558316 0 -10627.857 -10794.38 -10794.224 -393.99
4 14.050859 18.303344 0.27279397 0.95886179 -170.78053 -169.82167 3.1495694 0 -10627.837 -10794.509 -10794.236 -408.93357
5 19.424842 28.194049 0.42020553 0.9585778 -170.83354 -169.87496 3.016716 0 -10627.813 -10794.671 -10794.251 -428.10518
6 25.315002 39.895741 0.59460812 0.95823491 -170.90236 -169.94412 2.8599973 0 -10627.779 -10794.864 -10794.269 -451.5421
7 31.945826 53.187854 0.79271445 0.95783553 -170.98485 -170.02702 2.6826348 0 -10627.739 -10795.083 -10794.29 -479.20295
8 38.843432 67.819953 1.010792 0.95738252 -171.0742 -170.11682 2.4882829 0 -10627.697 -10795.325 -10794.314 -510.99723
9 46.718156 83.516497 1.2447341 0.95687915 -171.17882 -170.22194 2.2809581 0 -10627.645 -10795.586 -10794.342 -546.9588
10 54.366233 99.982124 1.4901386 0.95632904 -171.28928 -170.33295 2.0649603 0 -10627.594 -10795.862 -10794.372 -586.9531
Loop time of 54.3663 on 1 procs for 10 steps with 6 atoms
Performance: 0.004 ns/day, 6063.691 hours/ns, 0.458 timesteps/s, 2.749 atom-step/s
Performance: 0.002 ns/day, 15101.736 hours/ns, 0.184 timesteps/s, 1.104 atom-step/s
100.0% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.572e-06 | 7.572e-06 | 7.572e-06 | 0.0 | 0.00
Bond | 5.119e-06 | 5.119e-06 | 5.119e-06 | 0.0 | 0.00
Pair | 3.4059e-05 | 3.4059e-05 | 3.4059e-05 | 0.0 | 0.00
Bond | 1.1659e-05 | 1.1659e-05 | 1.1659e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 2.881e-06 | 2.881e-06 | 2.881e-06 | 0.0 | 0.00
Output | 5.439e-05 | 5.439e-05 | 5.439e-05 | 0.0 | 0.00
Modify | 4.3658 | 4.3658 | 4.3658 | 0.0 |100.00
Other | | 6.093e-06 | | | 0.00
Comm | 1.0583e-05 | 1.0583e-05 | 1.0583e-05 | 0.0 | 0.00
Output | 0.00016554 | 0.00016554 | 0.00016554 | 0.0 | 0.00
Modify | 54.366 | 54.366 | 54.366 | 0.0 |100.00
Other | | 1.567e-05 | | | 0.00
Nlocal: 6 ave 6 max 6 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -146,4 +154,4 @@ Ave neighs/atom = 2
Ave special neighs/atom = 1
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:01:11
Total wall time: 0:01:59

View File

@ -96,8 +96,8 @@ thermo 1
variable p equal extract_setting(world_size)
mdi plugin nwchem_mdi mdi "-role ENGINE -name NWChem -method LINK" extra "template.water.nw water.dimer.nw log.water.pwdft.qmmm.plugin.$p" command "run 1"
run 1
mdi plugin nwchem_mdi mdi "-role ENGINE -name NWChem -method LINK" extra "template.water.nw water.nw log.water.pwdft.qmmm.plugin.$p" command "run 10"
run 10
Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
@ -118,23 +118,32 @@ Neighbor list info ...
bin: standard
Per MPI rank memory allocation (min/avg/max) = 7.453 | 7.453 | 7.453 Mbytes
Step CPU Temp KinEng E_vdwl E_coul E_pair E_mol E_long f_2 PotEng TotEng Press
0 0 0 0 0.95937448 -170.68735 -169.72797 3.3962859 0 -10627.878 -10794.209 -10794.209 -374.81058
1 1.2206038 1.1713705 0.017458165 0.95934217 -170.68987 -169.73053 3.3804649 0 -10627.878 -10794.229 -10794.211 -376.91596
Loop time of 1.22062 on 2 procs for 1 steps with 6 atoms
0 0 0 0 0.95937448 -170.68735 -169.72798 3.3962859 0 -10627.878 -10794.209 -10794.209 -374.81058
1 1.2511997 1.1713705 0.017458165 0.95934217 -170.68987 -169.73053 3.3804649 0 -10627.878 -10794.229 -10794.211 -376.91596
2 3.8385452 4.6633575 0.069502914 0.95924532 -170.70698 -169.74773 3.3333254 0 -10627.871 -10794.286 -10794.216 -383.31463
3 7.7938711 10.410148 0.15515337 0.95908482 -170.73725 -169.77816 3.2558316 0 -10627.857 -10794.38 -10794.224 -393.99
4 12.421381 18.303344 0.27279396 0.95886179 -170.78051 -169.82164 3.1495694 0 -10627.837 -10794.509 -10794.236 -408.93334
5 17.269932 28.194049 0.42020553 0.9585778 -170.83354 -169.87496 3.016716 0 -10627.813 -10794.671 -10794.251 -428.10518
6 22.792492 39.895742 0.59460814 0.95823491 -170.90235 -169.94412 2.8599973 0 -10627.779 -10794.864 -10794.269 -451.54206
7 28.768848 53.187856 0.79271447 0.95783553 -170.98485 -170.02701 2.6826348 0 -10627.739 -10795.083 -10794.29 -479.20292
8 34.74838 67.819955 1.010792 0.95738253 -171.07427 -170.11688 2.4882829 0 -10627.697 -10795.325 -10794.314 -510.99784
9 41.409904 83.516502 1.2447341 0.95687915 -171.17899 -170.22211 2.2809581 0 -10627.645 -10795.586 -10794.342 -546.96024
10 48.305399 99.982132 1.4901387 0.95632904 -171.28935 -170.33302 2.0649603 0 -10627.594 -10795.862 -10794.372 -586.95371
Loop time of 48.3054 on 2 procs for 10 steps with 6 atoms
Performance: 0.007 ns/day, 3390.614 hours/ns, 0.819 timesteps/s, 4.916 atom-step/s
88.7% CPU use with 2 MPI tasks x no OpenMP threads
Performance: 0.002 ns/day, 13418.171 hours/ns, 0.207 timesteps/s, 1.242 atom-step/s
93.0% CPU use with 2 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.563e-06 | 4.9165e-06 | 5.27e-06 | 0.0 | 0.00
Bond | 1.201e-06 | 2.0745e-06 | 2.948e-06 | 0.0 | 0.00
Pair | 4.6638e-05 | 4.8108e-05 | 4.9577e-05 | 0.0 | 0.00
Bond | 1.5399e-05 | 2.157e-05 | 2.7741e-05 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 7.511e-06 | 8.0675e-06 | 8.624e-06 | 0.0 | 0.00
Output | 1.7987e-05 | 1.9748e-05 | 2.151e-05 | 0.0 | 0.00
Modify | 1.2206 | 1.2206 | 1.2206 | 0.0 |100.00
Other | | 7.148e-06 | | | 0.00
Comm | 7.9132e-05 | 8.3424e-05 | 8.7716e-05 | 0.0 | 0.00
Output | 0.0001854 | 0.00020373 | 0.00022205 | 0.0 | 0.00
Modify | 48.305 | 48.305 | 48.305 | 0.0 |100.00
Other | | 4.187e-05 | | | 0.00
Nlocal: 3 ave 4 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 1
@ -148,4 +157,4 @@ Ave neighs/atom = 2
Ave special neighs/atom = 1
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:01:00
Total wall time: 0:01:45

View File

@ -58,9 +58,9 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
addflag = 1;
every = 1;
connectflag = 1;
elements = nullptr;
mcflag = 0;
id_mcfix = nullptr;
elements = nullptr;
int iarg = 3;
while (iarg < narg) {
@ -96,12 +96,6 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else
error->all(FLERR, "Illegal fix mdi/qm command");
iarg += 2;
} else if (strcmp(arg[iarg],"mc") == 0) {
if (iarg+2 > narg) error->all(FLERR, "Illegal fix mdi/qm command");
mcflag = 1;
delete[] id_mcfix;
id_mcfix = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg], "elements") == 0) {
const char *symbols[] = {
@ -133,6 +127,13 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
}
iarg += ntypes + 1;
} else if (strcmp(arg[iarg],"mc") == 0) {
if (iarg+2 > narg) error->all(FLERR, "Illegal fix mdi/qm command");
mcflag = 1;
delete[] id_mcfix;
id_mcfix = utils::strdup(arg[iarg+1]);
iarg += 2;
} else
error->all(FLERR, "Illegal fix mdi/qm command");
}
@ -143,16 +144,16 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
global_freq = every;
extscalar = 1;
peratom_flag = 1;
size_peratom_cols = 3;
peratom_freq = every;
extvector = 0;
if (virialflag) {
vector_flag = 1;
size_vector = 6;
extvector = 0;
}
peratom_flag = 1;
size_peratom_cols = 3;
peratom_freq = every;
if (addflag) {
energy_global_flag = 1;
virial_global_flag = 1;
@ -202,13 +203,11 @@ FixMDIQM::FixMDIQM(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// initialize outputs
qm_energy = 0.0;
if (virialflag) {
for (int i = 0; i < 6; i++) {
qm_virial[i] = 0.0;
virial[i] = 0.0;
}
sumflag = 0;
}
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
@ -311,7 +310,7 @@ void FixMDIQM::init()
if (ierr) error->all(FLERR, "MDI: >NATOMS command check");
MPI_Bcast(&natoms_exists, 1, MPI_INT, 0, world);
ierr = MDI_Check_command_exists("@DEFAULT", ">NATOMS", mdicomm, &celldispl_exists);
ierr = MDI_Check_command_exists("@DEFAULT", ">CELL_DISPL", mdicomm, &celldispl_exists);
if (ierr) error->all(FLERR, "MDI: >CELL_DISPL command check");
MPI_Bcast(&celldispl_exists, 1, MPI_INT, 0, world);
@ -322,6 +321,10 @@ void FixMDIQM::init()
ierr = MDI_Check_command_exists("@DEFAULT", ">TYPES", mdicomm, &types_exists);
if (ierr) error->all(FLERR, "MDI: >TYPES command check");
MPI_Bcast(&types_exists, 1, MPI_INT, 0, world);
ierr = MDI_Check_command_exists("@DEFAULT", "<STRESS", mdicomm, &stress_exists);
if (ierr) error->all(FLERR, "MDI: <STRESS command check");
MPI_Bcast(&stress_exists, 1, MPI_INT, 0, world);
}
// extract pointers to MC variables from id_mcfix
@ -339,7 +342,7 @@ void FixMDIQM::init()
}
// determine whether a new vs incremental QM calc is needed
// new when first run or when
// new if first run or if
// atom count or elements/types or box has changed between runs
// otherwise incremental = subsequent run of same system
@ -385,6 +388,7 @@ void FixMDIQM::init()
if (eqm[i] != eqm_old[i]) new_system = 1;
memory->destroy(eqm_old);
}
} else if (types_exists) {
if (new_system) set_tqm();
else {
@ -558,7 +562,7 @@ void FixMDIQM::post_force(int vflag)
// qm_virial_symmetric = fix output for global QM virial
// note MDI defines virial tensor as intensive (divided by volume), LAMMPS does not
if (vflag && virialflag) {
if (vflag && virialflag && stress_exists) {
ierr = MDI_Send_command("<STRESS", mdicomm);
if (ierr) error->all(FLERR, "MDI: <STRESS command");
ierr = MDI_Recv(qm_virial, 9, MDI_DOUBLE, mdicomm);
@ -689,6 +693,11 @@ void FixMDIQM::reallocate()
int FixMDIQM::set_nqm()
{
// require 3*nqm be a small INT, so can MPI_Allreduce xqm
if (3*atom->natoms > MAXSMALLINT)
error->all(FLERR,"Fix mdi/qm has too many atoms");
int ncount = atom->natoms;
nexclude = 0;
@ -913,11 +922,9 @@ void FixMDIQM::send_natoms()
if (natoms_exists) {
ierr = MDI_Send_command(">NATOMS", mdicomm);
if (ierr) error->all(FLERR, "MDI: >NATOMS command");
int n = static_cast<int>(atom->natoms);
ierr = MDI_Send(&n, 1, MDI_INT, mdicomm);
ierr = MDI_Send(&nqm, 1, MDI_INT, mdicomm);
if (ierr) error->all(FLERR, "MDI: >NATOMS data");
} else {
ierr = MDI_Send_command("<NATOMS", mdicomm);
if (ierr) error->all(FLERR, "MDI: <NATOMS command");
@ -926,7 +933,7 @@ void FixMDIQM::send_natoms()
if (ierr) error->all(FLERR, "MDI: <NATOMS data");
MPI_Bcast(&n, 1, MPI_INT, 0, world);
if (n != atom->natoms)
if (n != nqm)
error->all(FLERR, "MDI: Engine has wrong atom count and does not support >NATOMS command");
}
}
@ -952,7 +959,7 @@ void FixMDIQM::send_elements()
int ierr = MDI_Send_command(">ELEMENTS", mdicomm);
if (ierr) error->all(FLERR, "MDI: >ELEMENTS command");
ierr = MDI_Send(eqm, nqm, MDI_INT, mdicomm);
if (ierr) error->all(FLERR, "MDI: >ELEMETNS data");
if (ierr) error->all(FLERR, "MDI: >ELEMENTS data");
}
/* ----------------------------------------------------------------------

View File

@ -62,6 +62,7 @@ class FixMDIQM : public Fix {
MDI_Comm mdicomm;
int natoms_exists,celldispl_exists,elements_exists,types_exists;
int stress_exists;
int nmax;

File diff suppressed because it is too large Load Diff

View File

@ -42,12 +42,14 @@ class FixMDIQMMM : public Fix {
void min_post_neighbor() override;
void min_pre_force(int) override;
void min_post_force(int) override;
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 compute_scalar() override;
double compute_vector(int) override;
double memory_usage() override;
private:
@ -59,49 +61,57 @@ class FixMDIQMMM : public Fix {
int *elements;
int mode; // QMMM method = DIRECT or POTENTIAL
int lmpunits; // REAL, METAL, or NATIVE
int first_send; // 1 until initial info passed to MDI engine
double qm_cell[9],qm_cell_displ[3];
double qm_energy;
double qm_virial[9], qm_virial_symmetric[6];
MDI_Comm mdicomm;
int natoms_exists,celldispl_exists,elements_exists,types_exists;
int stress_exists;
int nmax;
class Pair *pair_coul; // ptr to instance of pair coul variant
// data for QM portion
// QM atom data structs
int nqm; // # of QM atoms
int nqm_last,max_nqm;
tagint *qmIDs; // IDs of QM atoms in ascending order
double **xqm,**fqm; // QM coords and forces
double *qqm; // QM charges
int *eqm; // QM atom atomic numbers
double *qpotential; // Coulomb potential
double **xqm_mine; // same values for QM atoms I own
double *qqm_mine;
double *qpotential_mine;
int *eqm_mine;
int *qm2owned; // index of local atom for each QM atom
// index = -1 if this proc does not own
int *eqm,*eqm_mine;
int *tqm,*tqm_mine;
double **xqm,**xqm_mine;
double *qqm,*qqm_mine;
double *qpotential,*qpotential_mine;
double **fqm;
double *ecoul; // peratom Coulombic energy from LAMMPS
int ncoulmax; // length of ecoul
// data for MM portion
// MM atom data structs
int nmm; // # of MM atoms
int nmm_last,max_nmm;
tagint *mmIDs; // IDs of MM atoms in ascending order
double **xmm,**fmm; // MM coords and forces
double *qmm; // MM charges
int *emm; // MM atom atomic numbers
double **xmm_mine; // same values for MM atoms I own
double *qmm_mine;
int *emm_mine;
int *mm2owned; // index of local atom for each MM atom
// index = -1 if this proc does not own
int *emm,*emm_mine;
int *tmm,*tmm_mine;
double **xmm,**xmm_mine;
double *qmm,*qmm_mine;
double **fmm;
// unit conversion factors
int lmpunits; // REAL, METAL, or NATIVE
double lmp2mdi_length, mdi2lmp_length;
double lmp2mdi_energy, mdi2lmp_energy;
double lmp2mdi_force, mdi2lmp_force;
@ -112,22 +122,39 @@ class FixMDIQMMM : public Fix {
void post_force_direct(int);
void post_force_potential(int);
void reallocate_qm();
void reallocate_mm();
int set_nqm();
int set_nmm();
void create_qm_list();
void create_mm_list();
void set_qm2owned();
void set_mm2owned();
void set_box();
void set_xqm();
void set_eqm();
void set_tqm();
void set_qqm();
void set_xqm();
void set_emm();
void set_qmm();
void set_xmm();
void set_emm();
void set_tmm();
void set_qmm();
void send_natoms_qm();
void send_types_qm();
void send_elements_qm();
void send_box();
void send_natoms_mm();
void send_types_mm();
void send_elements_mm();
void send_charges_mm();
void unit_conversions();
};