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

ELECTRODE package

This commit is contained in:
Ludwig Ahrens 2022-02-09 12:11:22 +01:00
parent d51017ff50
commit efa5db4c58
134 changed files with 85002 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
.gitlab-ci.yml Normal file
View File

@ -0,0 +1,38 @@
stages:
- unit-tests
image: ubuntu
unit-tests:
stage: unit-tests
script:
- echo "Installing dependencies"
- apt-get update
- apt-get install -y cmake make gcc g++ gfortran python3 python3-numpy libblas-dev liblapack-dev
- base_dir=$PWD
- echo "Building lammps in directory $base_dir"
- mkdir build && cd build
- cmake ../cmake
- cmake -D ENABLE_TESTING=on .
- cmake -DCMAKE_BUILD_TYPE=Debug .
- cmake -D PKG_MOLECULE=yes .
- cmake -D PKG_KSPACE=yes .
- cmake -D PKG_RIGID=yes .
- cmake -D PKG_MISC=yes .
- cmake -D PKG_ELECTRODE=yes .
- cmake --build .
- cmake -D PKG_INTEL=yes .
- cmake --build .
# need to run unit tests w/o root privileges
- adduser tester
- chown -R tester $base_dir
- runuser -u tester ctest
- echo "running tests in examples/electrode subdirectory of $base_dir"
- lmpbin=$base_dir/build/lmp
- cd $base_dir/examples/electrode/check_amat/
- bash check.sh $lmpbin
- cd $base_dir/examples/electrode/check_etypes/
- bash check.sh $lmpbin
- cd $base_dir/examples/electrode/check_ffield/
- bash check.sh $lmpbin
- cd $base_dir/examples/electrode/check_intel/
- bash check.sh $lmpbin

8
cmake/CMakeLists.txt
View File

@ -184,6 +184,7 @@ set(STANDARD_PACKAGES
DPD-SMOOTH DPD-SMOOTH
DRUDE DRUDE
EFF EFF
ELECTRODE
EXTRA-COMPUTE EXTRA-COMPUTE
EXTRA-DUMP EXTRA-DUMP
EXTRA-FIX EXTRA-FIX
@ -344,6 +345,7 @@ pkg_depends(DIELECTRIC KSPACE)
pkg_depends(DIELECTRIC EXTRA-PAIR) pkg_depends(DIELECTRIC EXTRA-PAIR)
pkg_depends(CG-DNA MOLECULE) pkg_depends(CG-DNA MOLECULE)
pkg_depends(CG-DNA ASPHERE) pkg_depends(CG-DNA ASPHERE)
pkg_depends(ELECTRODE KSPACE)
# detect if we may enable OpenMP support by default # detect if we may enable OpenMP support by default
set(BUILD_OMP_DEFAULT OFF) set(BUILD_OMP_DEFAULT OFF)
@ -379,7 +381,7 @@ if(BUILD_OMP)
target_link_libraries(lammps PRIVATE OpenMP::OpenMP_CXX) target_link_libraries(lammps PRIVATE OpenMP::OpenMP_CXX)
endif() endif()
if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_LATTE) if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_LATTE OR PKG_ELECTRODE)
enable_language(C) enable_language(C)
find_package(LAPACK) find_package(LAPACK)
find_package(BLAS) find_package(BLAS)
@ -582,6 +584,10 @@ foreach(PKG_LIB POEMS ATC AWPMD H5MD MESONT)
endif() endif()
endforeach() endforeach()
if(PKG_ELECTRODE)
target_link_libraries(lammps PRIVATE ${LAPACK_LIBRARIES})
endif()
if(PKG_AWPMD) if(PKG_AWPMD)
target_link_libraries(awpmd PRIVATE ${LAPACK_LIBRARIES}) target_link_libraries(awpmd PRIVATE ${LAPACK_LIBRARIES})
endif() endif()

4
cmake/Modules/Packages/INTEL.cmake
View File

@ -112,5 +112,9 @@ if(PKG_KSPACE)
RegisterIntegrateStyle(${INTEL_SOURCES_DIR}/verlet_lrt_intel.h) RegisterIntegrateStyle(${INTEL_SOURCES_DIR}/verlet_lrt_intel.h)
endif() endif()
if(PKG_ELECTRODE)
list(APPEND INTEL_SOURCES ${INTEL_SOURCES_DIR}/electrode_accel_intel.cpp)
endif()
target_sources(lammps PRIVATE ${INTEL_SOURCES}) target_sources(lammps PRIVATE ${INTEL_SOURCES})
target_include_directories(lammps PRIVATE ${INTEL_SOURCES_DIR}) target_include_directories(lammps PRIVATE ${INTEL_SOURCES_DIR})

35
doc/src/Build_extras.rst
View File

@ -1289,6 +1289,41 @@ be built for the most part with all major versions of the C++ language.
---------- ----------
.. _electrode:
ELECTRODE package
-----------------
This package depends on the KSPACE package.
.. tabs::
.. tab:: CMake build
No additional settings are needed besides ``-D PKG_KSPACE=yes`` and ``-D
PKG_ELECTRODE=yes``.
.. tab:: Traditional make
The package is activated with ``make yes-KSPACE`` and ``make
yes-ELECTRODE``
Note that the ``Makefile.lammps`` file has settings for the BLAS and
LAPACK linear algebra libraries. As explained in ``lib/awpmd/README``
these can either exist on your system, or you can use the files provided
in ``lib/linalg``. In the latter case you also need to build the library
in ``lib/linalg`` with a command like these:
.. code-block:: bash
$ make lib-linalg # print help message
$ make lib-linalg args="-m serial" # build with GNU Fortran compiler (settings as with "make serial")
$ make lib-linalg args="-m mpi" # build with default MPI Fortran compiler (settings as with "make mpi")
$ make lib-linalg args="-m gfortran" # build with GNU Fortran compiler
----------
.. _ml-pace: .. _ml-pace:
ML-PACE package ML-PACE package

3
doc/src/Commands_fix.rst
View File

@ -66,6 +66,9 @@ OPT.
* :doc:`edpd/source <fix_dpd_source>` * :doc:`edpd/source <fix_dpd_source>`
* :doc:`efield <fix_efield>` * :doc:`efield <fix_efield>`
* :doc:`ehex <fix_ehex>` * :doc:`ehex <fix_ehex>`
* :doc:`electrode/conp <fix_electrode_conp>`
* :doc:`electrode/conq <fix_electrode_conq>`
* :doc:`electrode/thermo <fix_electrode_thermo>`
* :doc:`electron/stopping <fix_electron_stopping>` * :doc:`electron/stopping <fix_electron_stopping>`
* :doc:`electron/stopping/fit <fix_electron_stopping>` * :doc:`electron/stopping/fit <fix_electron_stopping>`
* :doc:`enforce2d (k) <fix_enforce2d>` * :doc:`enforce2d (k) <fix_enforce2d>`

2
doc/src/Commands_kspace.rst
View File

@ -27,6 +27,7 @@ OPT.
* :doc:`ewald/disp/dipole <kspace_style>` * :doc:`ewald/disp/dipole <kspace_style>`
* :doc:`ewald/dipole <kspace_style>` * :doc:`ewald/dipole <kspace_style>`
* :doc:`ewald/dipole/spin <kspace_style>` * :doc:`ewald/dipole/spin <kspace_style>`
* :doc:`ewald/electrode <kspace_style>`
* :doc:`msm (o) <kspace_style>` * :doc:`msm (o) <kspace_style>`
* :doc:`msm/cg (o) <kspace_style>` * :doc:`msm/cg (o) <kspace_style>`
* :doc:`msm/dielectric <kspace_style>` * :doc:`msm/dielectric <kspace_style>`
@ -41,4 +42,5 @@ OPT.
* :doc:`pppm/stagger <kspace_style>` * :doc:`pppm/stagger <kspace_style>`
* :doc:`pppm/tip4p (o) <kspace_style>` * :doc:`pppm/tip4p (o) <kspace_style>`
* :doc:`pppm/dielectric <kspace_style>` * :doc:`pppm/dielectric <kspace_style>`
* :doc:`pppm/electrode <kspace_style>`
* :doc:`scafacos <kspace_style>` * :doc:`scafacos <kspace_style>`

23
doc/src/Packages_details.rst
View File

@ -49,6 +49,7 @@ page gives those details.
* :ref:`DPD-SMOOTH <PKG-DPD-SMOOTH>` * :ref:`DPD-SMOOTH <PKG-DPD-SMOOTH>`
* :ref:`DRUDE <PKG-DRUDE>` * :ref:`DRUDE <PKG-DRUDE>`
* :ref:`EFF <PKG-EFF>` * :ref:`EFF <PKG-EFF>`
* :ref:`ELECTRODE <PKG-ELECTRODE>`
* :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>` * :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
* :ref:`EXTRA-DUMP <PKG-EXTRA-DUMP>` * :ref:`EXTRA-DUMP <PKG-EXTRA-DUMP>`
* :ref:`EXTRA-FIX <PKG-EXTRA-FIX>` * :ref:`EXTRA-FIX <PKG-EXTRA-FIX>`
@ -770,6 +771,28 @@ tools/eff; see its README file.
------------------- -------------------
.. _PKG-ELECTRODE:
ELECTRODE package
--------------------
**Contents:**
The ELECTRODE package allows the user to enforce a constant potential method for
groups of atoms that interact with the remaining atoms as electrolyte.
**Authors:** The ELECTRODE library is written and maintained by Ludwig Ahrens-Iwers
(TUHH, Hamburg, Germany), Shern Tee (UQ, Brisbane, Australia) and Robert Meissner (TUHH, Hamburg, Germany).
**Install:**
This package has :ref:`specific installation instructions <electrode>` on the :doc:`Build extras <Build_extras>` page.
**Supporting info:**
* :doc:`fix electrode/conp <fix_electrode_conp>`
----------
.. _PKG-EXTRA-COMPUTE: .. _PKG-EXTRA-COMPUTE:
EXTRA-COMPUTE package EXTRA-COMPUTE package

5
doc/src/Packages_list.rst
View File

@ -143,6 +143,11 @@ whether an extra library is needed to build and use the package:
- :doc:`pair_style eff/cut <pair_eff>` - :doc:`pair_style eff/cut <pair_eff>`
- PACKAGES/eff - PACKAGES/eff
- no - no
* - :ref:`ELECTRODE <PKG-ELECTRODE>`
- electrode charges to match potential
- :doc:`fix electrode/conp <fix_electrode_conp>`
- PACKAGES/electrode
- no
* - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>` * - :ref:`EXTRA-COMPUTE <PKG-EXTRA-COMPUTE>`
- additional compute styles - additional compute styles
- :doc:`compute <compute>` - :doc:`compute <compute>`

3
doc/src/fix.rst
View File

@ -209,6 +209,9 @@ accelerated styles exist.
* :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations * :doc:`edpd/source <fix_dpd_source>` - add heat source to eDPD simulations
* :doc:`efield <fix_efield>` - impose electric field on system * :doc:`efield <fix_efield>` - impose electric field on system
* :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm * :doc:`ehex <fix_ehex>` - enhanced heat exchange algorithm
* :doc:`electrode/conp <fix_electrode_conp>` - impose electric potential
* :doc:`electrode/conq <fix_electrode_conq>` - impose total electric charge
* :doc:`electrode/thermo <fix_electrode_thermo>` - apply thermo-potentiostat
* :doc:`electron/stopping <fix_electron_stopping>` - electronic stopping power as a friction force * :doc:`electron/stopping <fix_electron_stopping>` - electronic stopping power as a friction force
* :doc:`electron/stopping/fit <fix_electron_stopping>` - electronic stopping power as a friction force * :doc:`electron/stopping/fit <fix_electron_stopping>` - electronic stopping power as a friction force
* :doc:`enforce2d <fix_enforce2d>` - zero out z-dimension velocity and force * :doc:`enforce2d <fix_enforce2d>` - zero out z-dimension velocity and force

100
doc/src/fix_electrode_conp.rst Normal file
View File

@ -0,0 +1,100 @@
.. index:: fix electrode/conp
fix electrode/conp command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID charge_udpate potential eta keyword values ...
* ID = name of fix
* group-ID = name of group fix is applied to
* potential = electric potential in Volts
* eta = reciprocal width of electrode charge densities
.. parsed-literal::
*symm(etry) on/off*
turn on/off charge neutrality constraint
*couple <group-ID> <potential>*
group-ID = add another group of electrode atoms
potential = electric potential in volts applied to this electrode
*etypes <types>*
specify atom types exclusive to the electrode for optimized neighbor
lists
*ffield on/off*
turn on/off finite-field implementation
*write_mat <filename>*
write elastance matrix to file
*write_inv <filename>*
write inverted matrix to file
*read_mat <filename>*
read elastance matrix from file
*read_inv <filename>*
read inverted matrix from file
Examples
""""""""
.. code-block:: LAMMPS
fix fxconp electrodes electrode/conp 0.0 1.805 symm on
fix fxconp bot electrode/conp -1.0 1.805 couple top 1.0 write_inv inv.csv symm on
Description
"""""""""""
The fix electrode/conp implements a constant potential method (CPM)
(:ref:`Siepmann <Siepmann>`, :ref:`Reed <Reed3>`). Charges of groups specified
as group-ID and with the `couple` keyword are adapted to meet their respective
potential at every time step. An arbitrary number of electrodes can be set but
the respective groups may not overlap. Electrode charges have a Gaussian charge
distribution with reciprocal width eta. The energy minimization is achieved via
matrix inversion :ref:`(Wang) <Wang5>`.
The fix necessitates the use of a long range solver that can provide the matrix
of electrode-electrode interactions and a vector of electrode-electrolyte
interactions. The Kspace styles *ewald/electrode* and *pppm/electrode*
:ref:`(Ahrens-Iwers) <Ahrens-Iwers>` are created specifically for this task.
For systems with non-periodic boundaries in one or two directions dipole
corrections are available with the :doc:`kspace_modify <kspace_modify>`. For
ewald/electrode a two-dimensional Ewald summation :ref:`(Hu) <Hu>` can be used
by setting "slab ew2d":
.. code-block:: LAMMPS
kspace_modify slab <slab_factor>
kspace_modify wire <wire_factor>
kspace_modify slab ew2d
.. warning::
Atom positions of electrode particles have to be fixed at all times.
----------
.. _Siepmann:
**(Siepmann)** Siepmann and Strik, J. Chem. Phys. 102, 511 (1995).
.. _Reed3:
**(Reed)** Reed *et al.*, J. Chem. Phys. 126, 084704 (2007).
.. _Wang5:
**(Wang)** Wang *et al.*, J. Chem. Phys. 141, 184102 (2014).
.. _Ahrens-Iwers:
**(Ahrens-Iwers)** Ahrens-Iwers and Meissner, J. Chem. Phys. 155, 104104 (2021).
.. _Hu:
**(Hu)** Hu, J. Chem. Theory Comput. 10, 5254 (2014).

61
doc/src/fix_electrode_conq.rst Normal file
View File

@ -0,0 +1,61 @@
.. index:: fix electrode/conq
fix electrode/conq command
==========================
Syntax
""""""
.. parsed-literal::
fix ID group-ID charge_udpate potential eta keyword values ...
* ID = name of fix
* group-ID = name of group fix is applied to
* potential = electric potential in Volts
* eta = reciprocal width of electrode charge densities
.. parsed-literal::
*couple <group-ID> <potential>*
group-ID = add another group of electrode atoms
potential = electric potential in volts applied to this electrode
*etypes <types>*
specify atom types exclusive to the electrode for optimized neighbor
lists
*ffield on/off*
turn on/off finite-field implementation
*write_mat <filename>*
write elastance matrix to file
*write_inv <filename>*
write inverted matrix to file
*read_mat <filename>*
read elastance matrix from file
*read_inv <filename>*
read inverted matrix from file
Examples
""""""""
.. code-block:: LAMMPS
fix fxconq electrodes electrode/conq 0.0 1.805 symm on
fix fxconq bot electrode/conq -1.0 1.805 couple top 1.0 write_inv inv.csv symm on
Description
"""""""""""
Lorem ipsum
.. code-block:: LAMMPS
kspace_modify slab <slab_factor>
kspace_modify wire <wire_factor>
kspace_modify slab ew2d
.. warning::
Atom positions of electrode particles have to be fixed at all times.
----------

65
doc/src/fix_electrode_thermo.rst Normal file
View File

@ -0,0 +1,65 @@
.. index:: fix electrode/thermo
fix electrode/thermo command
============================
Syntax
""""""
.. parsed-literal::
fix ID group-ID charge_udpate potential eta temp T_v tau_v keyword values ...
* ID = name of fix
* group-ID = name of group fix is applied to
* potential = electric potential in Volts
* eta = reciprocal width of electrode charge densities
* T_v temperature parameter of thermo-potentiostat
* tau_v time parameter of thermo-potentiostat
.. parsed-literal::
*symm(etry) on/off*
turn on/off charge neutrality constraint
*couple <group-ID> <potential>*
group-ID = add another group of electrode atoms
potential = electric potential in volts applied to this electrode
*etypes <types>*
specify atom types exclusive to the electrode for optimized neighbor
lists
*ffield on/off*
turn on/off finite-field implementation
*write_mat <filename>*
write elastance matrix to file
*write_inv <filename>*
write inverted matrix to file
*read_mat <filename>*
read elastance matrix from file
*read_inv <filename>*
read inverted matrix from file
Examples
""""""""
.. code-block:: LAMMPS
fix fxthermo electrodes electrode/thermo 0.0 1.805 symm on
fix fxthermo bot electrode/thermo -1.0 1.805 couple top 1.0 write_inv inv.csv symm on
Description
"""""""""""
Lorem ipsum
.. code-block:: LAMMPS
kspace_modify slab <slab_factor>
kspace_modify wire <wire_factor>
kspace_modify slab ew2d
.. warning::
Atom positions of electrode particles have to be fixed at all times.
----------

8
doc/src/kspace_style.rst
View File

@ -4,6 +4,7 @@
.. index:: kspace_style ewald/disp .. index:: kspace_style ewald/disp
.. index:: kspace_style ewald/disp/dipole .. index:: kspace_style ewald/disp/dipole
.. index:: kspace_style ewald/omp .. index:: kspace_style ewald/omp
.. index:: kspace_style ewald/electrode
.. index:: kspace_style pppm .. index:: kspace_style pppm
.. index:: kspace_style pppm/kk .. index:: kspace_style pppm/kk
.. index:: kspace_style pppm/omp .. index:: kspace_style pppm/omp
@ -23,6 +24,7 @@
.. index:: kspace_style pppm/stagger .. index:: kspace_style pppm/stagger
.. index:: kspace_style pppm/tip4p .. index:: kspace_style pppm/tip4p
.. index:: kspace_style pppm/tip4p/omp .. index:: kspace_style pppm/tip4p/omp
.. index:: kspace_style pppm/electrode
.. index:: kspace_style msm .. index:: kspace_style msm
.. index:: kspace_style msm/omp .. index:: kspace_style msm/omp
.. index:: kspace_style msm/cg .. index:: kspace_style msm/cg
@ -40,7 +42,7 @@ Syntax
kspace_style style value kspace_style style value
* style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/disp/dipole* or *ewald/omp* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *pppm/disp/dielectric* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos* * style = *none* or *ewald* or *ewald/dipole* or *ewald/dipole/spin* or *ewald/disp* or *ewald/disp/dipole* or *ewald/omp* or *ewald/electrode* or *pppm* or *pppm/cg* or *pppm/disp* or *pppm/tip4p* or *pppm/stagger* or *pppm/disp/tip4p* or *pppm/gpu* or *pppm/intel* or *pppm/disp/intel* or *pppm/kk* or *pppm/omp* or *pppm/cg/omp* or *pppm/disp/tip4p/omp* or *pppm/tip4p/omp* or *pppm/dielectic* or *pppm/disp/dielectric* or *pppm/electrode* or *msm* or *msm/cg* or *msm/omp* or *msm/cg/omp* or *msm/dielectric* or *scafacos*
.. parsed-literal:: .. parsed-literal::
@ -57,6 +59,8 @@ Syntax
accuracy = desired relative error in forces accuracy = desired relative error in forces
*ewald/omp* value = accuracy *ewald/omp* value = accuracy
accuracy = desired relative error in forces accuracy = desired relative error in forces
*ewald/electrode* value = accuracy
accuracy = desired relative error in forces
*pppm* value = accuracy *pppm* value = accuracy
accuracy = desired relative error in forces accuracy = desired relative error in forces
*pppm/cg* values = accuracy (smallq) *pppm/cg* values = accuracy (smallq)
@ -97,6 +101,8 @@ Syntax
accuracy = desired relative error in forces accuracy = desired relative error in forces
*pppm/disp/dielectric* value = accuracy *pppm/disp/dielectric* value = accuracy
accuracy = desired relative error in forces accuracy = desired relative error in forces
*pppm/electrode* value = accuracy
accuracy = desired relative error in forces
*msm* value = accuracy *msm* value = accuracy
accuracy = desired relative error in forces accuracy = desired relative error in forces
*msm/cg* value = accuracy (smallq) *msm/cg* value = accuracy (smallq)

2
examples/electrode/check_amat/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
*.csv
log.lammps_*

600
examples/electrode/check_amat/charged.data Normal file
View File

@ -0,0 +1,600 @@
LAMMPS data file via write_data, version 24 Dec 2020, timestep = 0
288 atoms
2 atom types
0.285001 25.325067 xlo xhi
-0.368323 24.671743 ylo yhi
-1.043333 43.303354 zlo zhi
Masses
1 196.966553
2 196.966553
Pair Coeffs # lj/cut/coul/long
1 0 0
2 0 0
Atoms # full
2 1 1 -0.002921 2.086672 2.086672 8 0 0 0
3 1 1 -0.04033 2.086672 0 10.086673 0 1 0
5 1 1 -0.000757 4.173345 0 8 0 1 0
6 1 1 -0.000483 6.260016 2.086672 8 0 0 0
7 1 1 -0.02377 6.260016 0 10.086673 0 1 0
8 1 1 -0.031311 4.173345 2.086672 10.086673 0 0 0
9 1 1 -0.002302 8.346689 0 8 0 1 0
12 1 1 -0.03014 8.346689 2.086672 10.086673 0 0 0
26 1 1 -3e-05 2.086672 6.260016 8 0 0 0
27 1 1 -0.011104 2.086672 4.173345 10.086673 0 0 0
29 1 1 0.000807 4.173345 4.173345 8 0 0 0
30 1 1 -0.000974 6.260016 6.260016 8 0 0 0
31 1 1 -0.017105 6.260016 4.173345 10.086673 0 0 0
32 1 1 -0.010937 4.173345 6.260016 10.086673 0 0 0
33 1 1 0.001587 8.346689 4.173345 8 0 0 0
36 1 1 0.004269 8.346689 6.260016 10.086673 0 0 0
10 1 1 -0.000534 10.433361 2.086672 8 0 0 0
11 1 1 -0.006246 10.433361 0 10.086673 0 1 0
13 1 1 0.000222 12.520033 0 8 0 1 0
14 1 1 -0.000958 14.606705 2.086672 8 0 0 0
15 1 1 -0.024586 14.606705 0 10.086673 0 1 0
16 1 1 -0.032305 12.520033 2.086672 10.086673 0 0 0
17 1 1 -0.001102 16.693378 0 8 0 1 0
20 1 1 -0.015815 16.693378 2.086672 10.086673 0 0 0
34 1 1 0.001733 10.433361 6.260016 8 0 0 0
35 1 1 -0.001065 10.433361 4.173345 10.086673 0 0 0
37 1 1 -0.000623 12.520033 4.173345 8 0 0 0
38 1 1 -0.003876 14.606705 6.260016 8 0 0 0
39 1 1 -0.030292 14.606705 4.173345 10.086673 0 0 0
40 1 1 -0.024863 12.520033 6.260016 10.086673 0 0 0
41 1 1 -0.000446 16.693378 4.173345 8 0 0 0
44 1 1 -0.044616 16.693378 6.260016 10.086673 0 0 0
1 1 1 -0.00023 25.040066 0 8 0 1 0
4 1 1 -0.012754 25.040066 2.086672 10.086673 0 0 0
18 1 1 0.00046 18.780048 2.086672 8 0 0 0
19 1 1 -0.011993 18.780048 0 10.086673 0 1 0
21 1 1 0.001654 20.86672 0 8 0 1 0
22 1 1 0.000999 22.953394 2.086672 8 0 0 0
23 1 1 -0.010485 22.953394 0 10.086673 0 1 0
24 1 1 -0.038444 20.86672 2.086672 10.086673 0 0 0
25 1 1 3.8e-05 25.040066 4.173345 8 0 0 0
28 1 1 -0.032472 25.040066 6.260016 10.086673 0 0 0
42 1 1 -0.002036 18.780048 6.260016 8 0 0 0
43 1 1 -0.031052 18.780048 4.173345 10.086673 0 0 0
45 1 1 -0.001031 20.86672 4.173345 8 0 0 0
46 1 1 -0.00021 22.953394 6.260016 8 0 0 0
47 1 1 -0.004291 22.953394 4.173345 10.086673 0 0 0
48 1 1 -0.020295 20.86672 6.260016 10.086673 0 0 0
50 1 1 -0.00328 2.086672 10.433361 8 0 0 0
51 1 1 -0.037649 2.086672 8.346689 10.086673 0 0 0
53 1 1 -0.000151 4.173345 8.346689 8 0 0 0
54 1 1 -0.000508 6.260016 10.433361 8 0 0 0
55 1 1 -0.014809 6.260016 8.346689 10.086673 0 0 0
56 1 1 -0.016074 4.173345 10.433361 10.086673 0 0 0
57 1 1 0.00092 8.346689 8.346689 8 0 0 0
60 1 1 -0.037118 8.346689 10.433361 10.086673 0 0 0
74 1 1 -0.000342 2.086672 14.606705 8 0 0 0
75 1 1 -0.027442 2.086672 12.520033 10.086673 0 0 0
77 1 1 0.000536 4.173345 12.520033 8 0 0 0
78 1 1 -0.000368 6.260016 14.606705 8 0 0 0
79 1 1 -0.022015 6.260016 12.520033 10.086673 0 0 0
80 1 1 0.003048 4.173345 14.606705 10.086673 0 0 0
81 1 1 -0.002516 8.346689 12.520033 8 0 0 0
84 1 1 -0.004673 8.346689 14.606705 10.086673 0 0 0
58 1 1 -0.003585 10.433361 10.433361 8 0 0 0
59 1 1 -0.015447 10.433361 8.346689 10.086673 0 0 0
61 1 1 -0.000856 12.520033 8.346689 8 0 0 0
62 1 1 0.001309 14.606705 10.433361 8 0 0 0
63 1 1 -0.034936 14.606705 8.346689 10.086673 0 0 0
64 1 1 -0.022444 12.520033 10.433361 10.086673 0 0 0
65 1 1 -0.00144 16.693378 8.346689 8 0 0 0
68 1 1 -0.008916 16.693378 10.433361 10.086673 0 0 0
82 1 1 0.002034 10.433361 14.606705 8 0 0 0
83 1 1 -0.012341 10.433361 12.520033 10.086673 0 0 0
85 1 1 0.000168 12.520033 12.520033 8 0 0 0
86 1 1 0.000653 14.606705 14.606705 8 0 0 0
87 1 1 0.007996 14.606705 12.520033 10.086673 0 0 0
88 1 1 0.015522 12.520033 14.606705 10.086673 0 0 0
89 1 1 -2.2e-05 16.693378 12.520033 8 0 0 0
92 1 1 -0.001057 16.693378 14.606705 10.086673 0 0 0
49 1 1 -0.006345 25.040066 8.346689 8 0 0 0
52 1 1 -0.05504 25.040066 10.433361 10.086673 0 0 0
66 1 1 -0.004359 18.780048 10.433361 8 0 0 0
67 1 1 -0.02956 18.780048 8.346689 10.086673 0 0 0
69 1 1 0.0004 20.86672 8.346689 8 0 0 0
70 1 1 0.000209 22.953394 10.433361 8 0 0 0
71 1 1 -0.017199 22.953394 8.346689 10.086673 0 0 0
72 1 1 -0.03354 20.86672 10.433361 10.086673 0 0 0
73 1 1 -8e-05 25.040066 12.520033 8 0 0 0
76 1 1 -0.025027 25.040066 14.606705 10.086673 0 0 0
90 1 1 0.000716 18.780048 14.606705 8 0 0 0
91 1 1 -0.013279 18.780048 12.520033 10.086673 0 0 0
93 1 1 -0.000971 20.86672 12.520033 8 0 0 0
94 1 1 -0.000871 22.953394 14.606705 8 0 0 0
95 1 1 -0.014468 22.953394 12.520033 10.086673 0 0 0
96 1 1 -0.028593 20.86672 14.606705 10.086673 0 0 0
98 1 1 -0.001941 2.086672 18.780048 8 0 0 0
99 1 1 -0.015815 2.086672 16.693378 10.086673 0 0 0
101 1 1 0.000289 4.173345 16.693378 8 0 0 0
102 1 1 7.1e-05 6.260016 18.780048 8 0 0 0
103 1 1 -0.03149 6.260016 16.693378 10.086673 0 0 0
104 1 1 -0.033754 4.173345 18.780048 10.086673 0 0 0
105 1 1 -0.003025 8.346689 16.693378 8 0 0 0
108 1 1 -0.025007 8.346689 18.780048 10.086673 0 0 0
122 1 1 0.000623 2.086672 22.953394 8 0 0 0
123 1 1 -0.007516 2.086672 20.86672 10.086673 0 0 0
125 1 1 0.000878 4.173345 20.86672 8 0 0 0
126 1 1 -0.001576 6.260016 22.953394 8 0 0 0
127 1 1 -0.029629 6.260016 20.86672 10.086673 0 0 0
128 1 1 0.004728 4.173345 22.953394 10.086673 0 0 0
129 1 1 -0.003433 8.346689 20.86672 8 0 0 0
132 1 1 -0.043264 8.346689 22.953394 10.086673 0 0 0
106 1 1 -4.9e-05 10.433361 18.780048 8 0 0 0
107 1 1 -0.026304 10.433361 16.693378 10.086673 0 0 0
109 1 1 0.001185 12.520033 16.693378 8 0 0 0
110 1 1 -0.000778 14.606705 18.780048 8 0 0 0
111 1 1 -0.020266 14.606705 16.693378 10.086673 0 0 0
112 1 1 -0.010312 12.520033 18.780048 10.086673 0 0 0
113 1 1 -0.000613 16.693378 16.693378 8 0 0 0
116 1 1 -0.029068 16.693378 18.780048 10.086673 0 0 0
130 1 1 -0.000577 10.433361 22.953394 8 0 0 0
131 1 1 -0.026467 10.433361 20.86672 10.086673 0 0 0
133 1 1 -0.000921 12.520033 20.86672 8 0 0 0
134 1 1 -0.000361 14.606705 22.953394 8 0 0 0
135 1 1 -0.017852 14.606705 20.86672 10.086673 0 0 0
136 1 1 -0.021653 12.520033 22.953394 10.086673 0 0 0
137 1 1 0.000268 16.693378 20.86672 8 0 0 0
140 1 1 -0.01542 16.693378 22.953394 10.086673 0 0 0
97 1 1 -0.000311 25.040066 16.693378 8 0 0 0
100 1 1 -0.012408 25.040066 18.780048 10.086673 0 0 0
114 1 1 -0.000958 18.780048 18.780048 8 0 0 0
115 1 1 -0.036602 18.780048 16.693378 10.086673 0 0 0
117 1 1 -0.005554 20.86672 16.693378 8 0 0 0
118 1 1 0.000738 22.953394 18.780048 8 0 0 0
119 1 1 -0.027189 22.953394 16.693378 10.086673 0 0 0
120 1 1 -0.023111 20.86672 18.780048 10.086673 0 0 0
121 1 1 7.5e-05 25.040066 20.86672 8 0 0 0
124 1 1 -0.014682 25.040066 22.953394 10.086673 0 0 0
138 1 1 -0.000345 18.780048 22.953394 8 0 0 0
139 1 1 -0.00903 18.780048 20.86672 10.086673 0 0 0
141 1 1 0.000128 20.86672 20.86672 8 0 0 0
142 1 1 0.000334 22.953394 22.953394 8 0 0 0
143 1 1 -0.002285 22.953394 20.86672 10.086673 0 0 0
144 1 1 -0.00745 20.86672 22.953394 10.086673 0 0 0
146 2 2 0.022955 2.086672 2.086672 32.173344 0 0 0
147 2 2 0.000876 2.086672 0 34.260021 0 1 0
149 2 2 0.027435 4.173345 0 32.173344 0 1 0
150 2 2 0.020862 6.260016 2.086672 32.173344 0 0 0
151 2 2 0.001578 6.260016 0 34.260021 0 1 0
152 2 2 0.001014 4.173345 2.086672 34.260021 0 0 0
153 2 2 0.021009 8.346689 0 32.173344 0 1 0
156 2 2 0.001184 8.346689 2.086672 34.260021 0 0 0
170 2 2 0.008185 2.086672 6.260016 32.173344 0 0 0
171 2 2 0.001602 2.086672 4.173345 34.260021 0 0 0
173 2 2 0.022473 4.173345 4.173345 32.173344 0 0 0
174 2 2 0.03292 6.260016 6.260016 32.173344 0 0 0
175 2 2 0.002046 6.260016 4.173345 34.260021 0 0 0
176 2 2 0.000803 4.173345 6.260016 34.260021 0 0 0
177 2 2 0.032031 8.346689 4.173345 32.173344 0 0 0
180 2 2 0.002258 8.346689 6.260016 34.260021 0 0 0
154 2 2 0.01514 10.433361 2.086672 32.173344 0 0 0
155 2 2 0.001872 10.433361 0 34.260021 0 1 0
157 2 2 0.004752 12.520033 0 32.173344 0 1 0
158 2 2 0.005408 14.606705 2.086672 32.173344 0 0 0
159 2 2 0.000237 14.606705 0 34.260021 0 1 0
160 2 2 0.001509 12.520033 2.086672 34.260021 0 0 0
161 2 2 0.012702 16.693378 0 32.173344 0 1 0
164 2 2 0.000611 16.693378 2.086672 34.260021 0 0 0
178 2 2 0.02783 10.433361 6.260016 32.173344 0 0 0
179 2 2 0.001241 10.433361 4.173345 34.260021 0 0 0
181 2 2 0.019546 12.520033 4.173345 32.173344 0 0 0
182 2 2 0.016796 14.606705 6.260016 32.173344 0 0 0
183 2 2 0.00231 14.606705 4.173345 34.260021 0 0 0
184 2 2 0.001481 12.520033 6.260016 34.260021 0 0 0
185 2 2 0.014646 16.693378 4.173345 32.173344 0 0 0
188 2 2 0.000878 16.693378 6.260016 34.260021 0 0 0
145 2 2 0.019659 25.040066 0 32.173344 0 1 0
148 2 2 0.00156 25.040066 2.086672 34.260021 0 0 0
162 2 2 0.017746 18.780048 2.086672 32.173344 0 0 0
163 2 2 0.001489 18.780048 0 34.260021 0 1 0
165 2 2 0.034109 20.86672 0 32.173344 0 1 0
166 2 2 0.031678 22.953394 2.086672 32.173344 0 0 0
167 2 2 0.002473 22.953394 0 34.260021 0 1 0
168 2 2 0.002352 20.86672 2.086672 34.260021 0 0 0
169 2 2 0.014267 25.040066 4.173345 32.173344 0 0 0
172 2 2 0.000168 25.040066 6.260016 34.260021 0 0 0
186 2 2 0.008992 18.780048 6.260016 32.173344 0 0 0
187 2 2 0.000504 18.780048 4.173345 34.260021 0 0 0
189 2 2 0.021772 20.86672 4.173345 32.173344 0 0 0
190 2 2 0.01304 22.953394 6.260016 32.173344 0 0 0
191 2 2 0.00118 22.953394 4.173345 34.260021 0 0 0
192 2 2 0.001457 20.86672 6.260016 34.260021 0 0 0
194 2 2 0.025689 2.086672 10.433361 32.173344 0 0 0
195 2 2 0.000649 2.086672 8.346689 34.260021 0 0 0
197 2 2 0.015557 4.173345 8.346689 32.173344 0 0 0
198 2 2 0.025584 6.260016 10.433361 32.173344 0 0 0
199 2 2 0.002292 6.260016 8.346689 34.260021 0 0 0
200 2 2 0.000819 4.173345 10.433361 34.260021 0 0 0
201 2 2 0.035591 8.346689 8.346689 32.173344 0 0 0
204 2 2 0.002678 8.346689 10.433361 34.260021 0 0 0
218 2 2 0.021343 2.086672 14.606705 32.173344 0 0 0
219 2 2 0.001799 2.086672 12.520033 34.260021 0 0 0
221 2 2 0.022952 4.173345 12.520033 32.173344 0 0 0
222 2 2 0.012325 6.260016 14.606705 32.173344 0 0 0
223 2 2 0.001876 6.260016 12.520033 34.260021 0 0 0
224 2 2 0.001718 4.173345 14.606705 34.260021 0 0 0
225 2 2 0.023719 8.346689 12.520033 32.173344 0 0 0
228 2 2 0.000744 8.346689 14.606705 34.260021 0 0 0
202 2 2 0.015873 10.433361 10.433361 32.173344 0 0 0
203 2 2 0.000758 10.433361 8.346689 34.260021 0 0 0
205 2 2 0.008097 12.520033 8.346689 32.173344 0 0 0
206 2 2 0.015646 14.606705 10.433361 32.173344 0 0 0
207 2 2 0.001702 14.606705 8.346689 34.260021 0 0 0
208 2 2 0.000994 12.520033 10.433361 34.260021 0 0 0
209 2 2 0.017819 16.693378 8.346689 32.173344 0 0 0
212 2 2 0.001677 16.693378 10.433361 34.260021 0 0 0
226 2 2 0.004714 10.433361 14.606705 32.173344 0 0 0
227 2 2 0.001251 10.433361 12.520033 34.260021 0 0 0
229 2 2 0.018756 12.520033 12.520033 32.173344 0 0 0
230 2 2 -0.00689 14.606705 14.606705 32.173344 0 0 0
231 2 2 0.000922 14.606705 12.520033 34.260021 0 0 0
232 2 2 -0.002432 12.520033 14.606705 34.260021 0 0 0
233 2 2 0.015625 16.693378 12.520033 32.173344 0 0 0
236 2 2 0.001044 16.693378 14.606705 34.260021 0 0 0
193 2 2 0.001698 25.040066 8.346689 32.173344 0 0 0
196 2 2 0.00215 25.040066 10.433361 34.260021 0 0 0
210 2 2 0.029386 18.780048 10.433361 32.173344 0 0 0
211 2 2 0.000676 18.780048 8.346689 34.260021 0 0 0
213 2 2 0.020182 20.86672 8.346689 32.173344 0 0 0
214 2 2 0.02881 22.953394 10.433361 32.173344 0 0 0
215 2 2 0.002276 22.953394 8.346689 34.260021 0 0 0
216 2 2 0.003118 20.86672 10.433361 34.260021 0 0 0
217 2 2 0.012257 25.040066 12.520033 32.173344 0 0 0
220 2 2 0.001591 25.040066 14.606705 34.260021 0 0 0
234 2 2 -0.009585 18.780048 14.606705 32.173344 0 0 0
235 2 2 0.001459 18.780048 12.520033 34.260021 0 0 0
237 2 2 0.020359 20.86672 12.520033 32.173344 0 0 0
238 2 2 0.004616 22.953394 14.606705 32.173344 0 0 0
239 2 2 0.000697 22.953394 12.520033 34.260021 0 0 0
240 2 2 5e-06 20.86672 14.606705 34.260021 0 0 0
242 2 2 0.034403 2.086672 18.780048 32.173344 0 0 0
243 2 2 0.002402 2.086672 16.693378 34.260021 0 0 0
245 2 2 0.021827 4.173345 16.693378 32.173344 0 0 0
246 2 2 0.011918 6.260016 18.780048 32.173344 0 0 0
247 2 2 0.000792 6.260016 16.693378 34.260021 0 0 0
248 2 2 0.002071 4.173345 18.780048 34.260021 0 0 0
249 2 2 -0.002694 8.346689 16.693378 32.173344 0 0 0
252 2 2 0.000351 8.346689 18.780048 34.260021 0 0 0
266 2 2 0.037372 2.086672 22.953394 32.173344 0 0 0
267 2 2 0.002492 2.086672 20.86672 34.260021 0 0 0
269 2 2 0.03734 4.173345 20.86672 32.173344 0 0 0
270 2 2 0.015818 6.260016 22.953394 32.173344 0 0 0
271 2 2 -0.0001 6.260016 20.86672 34.260021 0 0 0
272 2 2 0.002801 4.173345 22.953394 34.260021 0 0 0
273 2 2 0.017084 8.346689 20.86672 32.173344 0 0 0
276 2 2 0.001644 8.346689 22.953394 34.260021 0 0 0
250 2 2 -0.000573 10.433361 18.780048 32.173344 0 0 0
251 2 2 0.000448 10.433361 16.693378 34.260021 0 0 0
253 2 2 0.001141 12.520033 16.693378 32.173344 0 0 0
254 2 2 0.015575 14.606705 18.780048 32.173344 0 0 0
255 2 2 0.000258 14.606705 16.693378 34.260021 0 0 0
256 2 2 0.000372 12.520033 18.780048 34.260021 0 0 0
257 2 2 0.011285 16.693378 16.693378 32.173344 0 0 0
260 2 2 0.001199 16.693378 18.780048 34.260021 0 0 0
274 2 2 0.031072 10.433361 22.953394 32.173344 0 0 0
275 2 2 -0.00031 10.433361 20.86672 34.260021 0 0 0
277 2 2 -0.004414 12.520033 20.86672 32.173344 0 0 0
278 2 2 0.017397 14.606705 22.953394 32.173344 0 0 0
279 2 2 0.001344 14.606705 20.86672 34.260021 0 0 0
280 2 2 0.000672 12.520033 22.953394 34.260021 0 0 0
281 2 2 0.0252 16.693378 20.86672 32.173344 0 0 0
284 2 2 0.001824 16.693378 22.953394 34.260021 0 0 0
241 2 2 0.029057 25.040066 16.693378 32.173344 0 0 0
244 2 2 0.001076 25.040066 18.780048 34.260021 0 0 0
258 2 2 0.035183 18.780048 18.780048 32.173344 0 0 0
259 2 2 0.0002 18.780048 16.693378 34.260021 0 0 0
261 2 2 0.034083 20.86672 16.693378 32.173344 0 0 0
262 2 2 0.036995 22.953394 18.780048 32.173344 0 0 0
263 2 2 0.001484 22.953394 16.693378 34.260021 0 0 0
264 2 2 0.00377 20.86672 18.780048 34.260021 0 0 0
265 2 2 0.024028 25.040066 20.86672 32.173344 0 0 0
268 2 2 0.001334 25.040066 22.953394 34.260021 0 0 0
282 2 2 0.030927 18.780048 22.953394 32.173344 0 0 0
283 2 2 0.002185 18.780048 20.86672 34.260021 0 0 0
285 2 2 0.043618 20.86672 20.86672 32.173344 0 0 0
286 2 2 0.030615 22.953394 22.953394 32.173344 0 0 0
287 2 2 0.001068 22.953394 20.86672 34.260021 0 0 0
288 2 2 0.001887 20.86672 22.953394 34.260021 0 0 0
Velocities
2 0 0 0
3 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
12 0 0 0
26 0 0 0
27 0 0 0
29 0 0 0
30 0 0 0
31 0 0 0
32 0 0 0
33 0 0 0
36 0 0 0
10 0 0 0
11 0 0 0
13 0 0 0
14 0 0 0
15 0 0 0
16 0 0 0
17 0 0 0
20 0 0 0
34 0 0 0
35 0 0 0
37 0 0 0
38 0 0 0
39 0 0 0
40 0 0 0
41 0 0 0
44 0 0 0
1 0 0 0
4 0 0 0
18 0 0 0
19 0 0 0
21 0 0 0
22 0 0 0
23 0 0 0
24 0 0 0
25 0 0 0
28 0 0 0
42 0 0 0
43 0 0 0
45 0 0 0
46 0 0 0
47 0 0 0
48 0 0 0
50 0 0 0
51 0 0 0
53 0 0 0
54 0 0 0
55 0 0 0
56 0 0 0
57 0 0 0
60 0 0 0
74 0 0 0
75 0 0 0
77 0 0 0
78 0 0 0
79 0 0 0
80 0 0 0
81 0 0 0
84 0 0 0
58 0 0 0
59 0 0 0
61 0 0 0
62 0 0 0
63 0 0 0
64 0 0 0
65 0 0 0
68 0 0 0
82 0 0 0
83 0 0 0
85 0 0 0
86 0 0 0
87 0 0 0
88 0 0 0
89 0 0 0
92 0 0 0
49 0 0 0
52 0 0 0
66 0 0 0
67 0 0 0
69 0 0 0
70 0 0 0
71 0 0 0
72 0 0 0
73 0 0 0
76 0 0 0
90 0 0 0
91 0 0 0
93 0 0 0
94 0 0 0
95 0 0 0
96 0 0 0
98 0 0 0
99 0 0 0
101 0 0 0
102 0 0 0
103 0 0 0
104 0 0 0
105 0 0 0
108 0 0 0
122 0 0 0
123 0 0 0
125 0 0 0
126 0 0 0
127 0 0 0
128 0 0 0
129 0 0 0
132 0 0 0
106 0 0 0
107 0 0 0
109 0 0 0
110 0 0 0
111 0 0 0
112 0 0 0
113 0 0 0
116 0 0 0
130 0 0 0
131 0 0 0
133 0 0 0
134 0 0 0
135 0 0 0
136 0 0 0
137 0 0 0
140 0 0 0
97 0 0 0
100 0 0 0
114 0 0 0
115 0 0 0
117 0 0 0
118 0 0 0
119 0 0 0
120 0 0 0
121 0 0 0
124 0 0 0
138 0 0 0
139 0 0 0
141 0 0 0
142 0 0 0
143 0 0 0
144 0 0 0
146 0 0 0
147 0 0 0
149 0 0 0
150 0 0 0
151 0 0 0
152 0 0 0
153 0 0 0
156 0 0 0
170 0 0 0
171 0 0 0
173 0 0 0
174 0 0 0
175 0 0 0
176 0 0 0
177 0 0 0
180 0 0 0
154 0 0 0
155 0 0 0
157 0 0 0
158 0 0 0
159 0 0 0
160 0 0 0
161 0 0 0
164 0 0 0
178 0 0 0
179 0 0 0
181 0 0 0
182 0 0 0
183 0 0 0
184 0 0 0
185 0 0 0
188 0 0 0
145 0 0 0
148 0 0 0
162 0 0 0
163 0 0 0
165 0 0 0
166 0 0 0
167 0 0 0
168 0 0 0
169 0 0 0
172 0 0 0
186 0 0 0
187 0 0 0
189 0 0 0
190 0 0 0
191 0 0 0
192 0 0 0
194 0 0 0
195 0 0 0
197 0 0 0
198 0 0 0
199 0 0 0
200 0 0 0
201 0 0 0
204 0 0 0
218 0 0 0
219 0 0 0
221 0 0 0
222 0 0 0
223 0 0 0
224 0 0 0
225 0 0 0
228 0 0 0
202 0 0 0
203 0 0 0
205 0 0 0
206 0 0 0
207 0 0 0
208 0 0 0
209 0 0 0
212 0 0 0
226 0 0 0
227 0 0 0
229 0 0 0
230 0 0 0
231 0 0 0
232 0 0 0
233 0 0 0
236 0 0 0
193 0 0 0
196 0 0 0
210 0 0 0
211 0 0 0
213 0 0 0
214 0 0 0
215 0 0 0
216 0 0 0
217 0 0 0
220 0 0 0
234 0 0 0
235 0 0 0
237 0 0 0
238 0 0 0
239 0 0 0
240 0 0 0
242 0 0 0
243 0 0 0
245 0 0 0
246 0 0 0
247 0 0 0
248 0 0 0
249 0 0 0
252 0 0 0
266 0 0 0
267 0 0 0
269 0 0 0
270 0 0 0
271 0 0 0
272 0 0 0
273 0 0 0
276 0 0 0
250 0 0 0
251 0 0 0
253 0 0 0
254 0 0 0
255 0 0 0
256 0 0 0
257 0 0 0
260 0 0 0
274 0 0 0
275 0 0 0
277 0 0 0
278 0 0 0
279 0 0 0
280 0 0 0
281 0 0 0
284 0 0 0
241 0 0 0
244 0 0 0
258 0 0 0
259 0 0 0
261 0 0 0
262 0 0 0
263 0 0 0
264 0 0 0
265 0 0 0
268 0 0 0
282 0 0 0
283 0 0 0
285 0 0 0
286 0 0 0
287 0 0 0
288 0 0 0

17
examples/electrode/check_amat/charged.settings Normal file
View File

@ -0,0 +1,17 @@
atom_style full
units real
kspace_style pppm/electrode 1.0e-8
kspace_modify gewald 0.231
kspace_modify mesh 8 8 16
kspace_modify order 5
boundary p p f
kspace_modify slab 3.0
pair_style lj/cut/coul/long 12 12
processors * * 1
read_data "charged.data"
group ele type 1 2
pair_coeff * * 0.0 0.0

12
examples/electrode/check_amat/check.py Executable file
View File

@ -0,0 +1,12 @@
#!/bin/python3
import numpy as np
onestep = np.loadtxt("onestep.csv", skiprows=1)
twostep = np.loadtxt("twostep.csv", skiprows=1)
for matrix in [onestep, twostep]:
assert matrix.shape == (288, 288), matrix.shape
diff = abs(np.sum(onestep - twostep))
assert diff < 1e-11, diff

20
examples/electrode/check_amat/check.sh Executable file
View File

@ -0,0 +1,20 @@
#!/bin/bash -e
lmpbin=$1
if [ ! -f $lmpbin ]; then
echo "LAMMPS binary '$lmpbin' is not a file"
exit 1
fi
rm -f *.csv
$lmpbin < onestep.in > /dev/null
$lmpbin < twostep.in > /dev/null
for matrix in twostep.csv onestep.csv; do
if [ ! -f $matrix ]; then
echo "$matrix missing after running LAMMPS"
exit 1
fi
done
python3 check.py
rm *.csv log.lammps

6
examples/electrode/check_amat/onestep.in Normal file
View File

@ -0,0 +1,6 @@
include charged.settings
kspace_modify amat onestep
fix fxupdate ele electrode/conp 0.0 1.805 write_mat onestep.csv
run 0

6
examples/electrode/check_amat/twostep.in Normal file
View File

@ -0,0 +1,6 @@
include charged.settings
kspace_modify amat twostep
fix fxupdate ele electrode/conp 0.0 1.805 write_mat twostep.csv
run 0

2
examples/electrode/check_conq/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
log*
*csv

4766
examples/electrode/check_conq/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

19
examples/electrode/check_conq/cgil.settings Normal file
View File

@ -0,0 +1,19 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4

35
examples/electrode/check_conq/equal.in Normal file
View File

@ -0,0 +1,35 @@
variable tag string equal
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable qtop equal 1.0
variable qbot equal -v_qtop
variable vbot internal 0.0
variable vtop internal 0.0
fix fxupdate bot electrode/conq v_qbot 1.979 couple top v_qtop symm on
fix_modify fxupdate set v bot vbot
fix_modify fxupdate set v top vtop
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e) $(v_vbot:%20.15e) $(v_vtop:%20.15e)" file ${tag}.csv screen no title "#Fix electrode conq test setting = equal"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall v_vbot v_vtop
run 50

39
examples/electrode/check_conq/equal_intvars.in Normal file
View File

@ -0,0 +1,39 @@
variable tag string equal_intvars
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable qtop equal 1.0
variable qbot equal -v_qtop
variable qsb_bot internal 0.0
variable me00 internal 0.0
variable vbot equal v_me00*(v_qbot-v_qsb_bot)
variable vtop equal -v_vbot
fix fxupdate bot electrode/conp v_vbot 1.979 couple top v_vtop symm on
fix_modify fxupdate set qsb bot qsb_bot
fix_modify fxupdate set me bot bot me00
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e) $(v_vtop:%20.15e) $(v_vbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode conq test setting = equal, emulated with conp"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall v_vbot v_vtop
run 50

5
examples/electrode/check_conq/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0
#kspace_modify mesh 24 24 210
#kspace_modify order 5

31
examples/electrode/check_conq/ramp.in Normal file
View File

@ -0,0 +1,31 @@
variable tag string equal
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable vtop equal ramp(0,5)
variable vbot equal -v_vtop
fix fxupdate bot electrode/conq v_vbot 1.979 couple top v_vtop symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode conq test setting = ramp"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

2
examples/electrode/check_equalstyle/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
*.csv
log.lammps_*

4766
examples/electrode/check_equalstyle/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

19
examples/electrode/check_equalstyle/cgil.settings Normal file
View File

@ -0,0 +1,19 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4

20
examples/electrode/check_equalstyle/check.py Executable file
View File

@ -0,0 +1,20 @@
#!/bin/python3
import numpy as np
const = np.loadtxt("const.csv", skiprows=1)
equal = np.loadtxt("equal.csv", skiprows=1)
ramp = np.loadtxt("ramp.csv", skiprows=1)
for matrix in [const, equal, ramp]:
assert matrix.shape == (6, 3), matrix.shape
# equal follows const
diff0 = abs(np.sum(const - equal))
assert diff0 < 1e-8, diff0
# equal starts same as ramp but ends differently
diff1 = abs(np.sum(equal[0] - ramp[0]))
assert diff1 < 1e-8, diff1
diff2 = np.sum(abs(equal[-1] - ramp[-1]))
assert diff2 > 1e-8, diff2

22
examples/electrode/check_equalstyle/check.sh Executable file
View File

@ -0,0 +1,22 @@
#!/bin/bash -e
lmpbin=$1
if [ ! -f $lmpbin ]; then
echo "LAMMPS binary '$lmpbin' is not a file"
exit 1
fi
rm -f *.csv
$lmpbin < const.in > /dev/null
$lmpbin < equal.in > /dev/null
$lmpbin < ramp.in > /dev/null
for printout in const.csv equal.csv ramp.csv; do
if [ ! -f $printout ]; then
echo "$printout missing after running LAMMPS"
exit 1
fi
cat $printout
done
python3 check.py
rm *.csv log.lammps*

29
examples/electrode/check_equalstyle/const.in Normal file
View File

@ -0,0 +1,29 @@
variable tag string const
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = const"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

31
examples/electrode/check_equalstyle/equal.in Normal file
View File

@ -0,0 +1,31 @@
variable tag string equal
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable vtop equal 1.0
variable vbot equal -v_vtop
fix fxupdate bot electrode/conp v_vbot 1.979 couple top v_vtop symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = equal"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

5
examples/electrode/check_equalstyle/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0
#kspace_modify mesh 24 24 210
#kspace_modify order 5

31
examples/electrode/check_equalstyle/ramp.in Normal file
View File

@ -0,0 +1,31 @@
variable tag string ramp
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable vtop equal ramp(1.0,5.0)
variable vbot equal -v_vtop
fix fxupdate bot electrode/conp v_vbot 1.979 couple top v_vtop symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ramp"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

2
examples/electrode/check_etypes/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
log*
*csv

4766
examples/electrode/check_etypes/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

19
examples/electrode/check_etypes/cgil.settings Normal file
View File

@ -0,0 +1,19 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4

12
examples/electrode/check_etypes/check.py Executable file
View File

@ -0,0 +1,12 @@
#!/bin/python3
import numpy as np
ref = np.loadtxt("ref.csv", skiprows=1)
etypes = np.loadtxt("etypes.csv", skiprows=1)
for matrix in [ref, etypes]:
assert matrix.shape == (6, 3), matrix.shape
diff0 = abs(np.sum(ref - etypes))
assert diff0 < 1e-8, diff0

21
examples/electrode/check_etypes/check.sh Executable file
View File

@ -0,0 +1,21 @@
#!/bin/bash -e
lmpbin=$1
if [ ! -f $lmpbin ]; then
echo "LAMMPS binary '$lmpbin' is not a file"
exit 1
fi
rm -f *.csv
$lmpbin -i ref.in > /dev/null
$lmpbin -i etypes.in > /dev/null
for printout in ref.csv etypes.csv; do
if [ ! -f $printout ]; then
echo "$printout missing after running LAMMPS"
exit 1
fi
#cat $printout
done
python3 check.py
rm *.csv log.lammps*

31
examples/electrode/check_etypes/etypes.in Normal file
View File

@ -0,0 +1,31 @@
variable tag string etypes
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
neigh_modify check no # force neighbor list rebuild for cost testing
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on etypes 5
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = etypes"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

5
examples/electrode/check_etypes/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0 amat onestep # default
#kspace_modify mesh 24 24 210
#kspace_modify order 5

31
examples/electrode/check_etypes/ref.in Normal file
View File

@ -0,0 +1,31 @@
variable tag string ref
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
neigh_modify check no # force neighbor list rebuild for cost testing
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ref"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

2
examples/electrode/check_ffield/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
log*
*csv

4766
examples/electrode/check_ffield/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

21
examples/electrode/check_ffield/cgil.settings Normal file
View File

@ -0,0 +1,21 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4
neigh_modify one 4000

18
examples/electrode/check_ffield/check.py Executable file
View File

@ -0,0 +1,18 @@
#!/bin/python3
import numpy as np
ref = np.loadtxt("ref.csv", skiprows=1)
ffield = np.loadtxt("ffield.csv", skiprows=1)
ffield_flip = np.loadtxt("ffield_flip.csv", skiprows=1)
for matrix in [ref, ffield, ffield_flip]:
assert matrix.shape == (6, 3), matrix.shape
# ref is close to ffield
diff = np.sum(np.abs(ref - ffield))
assert diff < 3e-4, diff
# ffield and ffield_flip are identical
diff = np.sum(np.abs(ffield_flip - ffield))
assert diff < 1e-9, diff

22
examples/electrode/check_ffield/check.sh Executable file
View File

@ -0,0 +1,22 @@
#!/bin/bash -e
lmpbin=$1
if [ ! -f $lmpbin ]; then
echo "LAMMPS binary '$lmpbin' is not a file"
exit 1
fi
rm -f *.csv
$lmpbin < ref.in > /dev/null
$lmpbin < ffield.in > /dev/null
$lmpbin < ffield_flip.in > /dev/null
for printout in ref.csv ffield.csv ffield_flip.csv; do
if [ ! -f $printout ]; then
echo "$printout missing after running LAMMPS"
exit 1
fi
# cat $printout
done
python3 check.py
rm *.csv log.lammps*

29
examples/electrode/check_ffield/ffield.in Normal file
View File

@ -0,0 +1,29 @@
variable tag string ffield
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m.settings"
boundary p p p
pair_style lj/cut/coul/long 14 14
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on write_inv ${tag}_inv.csv ffield on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ffield"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

29
examples/electrode/check_ffield/ffield_flip.in Normal file
View File

@ -0,0 +1,29 @@
variable tag string ffield_flip
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m.settings"
boundary p p p
pair_style lj/cut/coul/long 14 14
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate top electrode/conp 1.0 1.979 couple bot -1.0 symm on write_inv ${tag}_inv.csv ffield on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ffield"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

4
examples/electrode/check_ffield/p3m.settings Normal file
View File

@ -0,0 +1,4 @@
kspace_style pppm/electrode 1.0e-6
#kspace_modify gewald 0.2311815
#kspace_modify mesh 24 24 210
#kspace_modify order 5

5
examples/electrode/check_ffield/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-6
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0
#kspace_modify mesh 24 24 210
#kspace_modify order 5

29
examples/electrode/check_ffield/ref.in Normal file
View File

@ -0,0 +1,29 @@
variable tag string ref
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 14 14
#processors * * 2
read_data "cap.data"
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on write_inv ${tag}_inv.csv
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ref"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

2
examples/electrode/check_hybrid/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
*.csv
log.lammps*

4766
examples/electrode/check_hybrid/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

17
examples/electrode/check_hybrid/cgil.settings Normal file
View File

@ -0,0 +1,17 @@
#Roy-Maroncelli BMI-PF6 + carbon
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4
set group bot charge 0
set group top charge 0
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on

32
examples/electrode/check_hybrid/coul.in Normal file
View File

@ -0,0 +1,32 @@
variable tag string coul
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style coul/long 16
#processors * * 2
read_data "cap.data"
pair_modify mix arithmetic
pair_coeff * *
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = coul"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

37
examples/electrode/check_hybrid/hybrid.in Normal file
View File

@ -0,0 +1,37 @@
variable tag string hybrid
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style hybrid/overlay lj/cut 16 coul/long 16
#processors * * 2
read_data "cap.data"
pair_modify mix arithmetic
pair_coeff 1 1 lj/cut 0.61185 4.38
pair_coeff 2 2 lj/cut 0.08604 3.41
pair_coeff 3 3 lj/cut 0.43738 5.04
pair_coeff 4 4 lj/cut 1.12572 5.06
pair_coeff 5 5 lj/cut 0.05497 3.37
pair_coeff * * coul/long
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = hybrid"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

36
examples/electrode/check_hybrid/lj.in Normal file
View File

@ -0,0 +1,36 @@
variable tag string lj
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
#include "p3m3dc.settings"
boundary p p f
pair_style lj/cut 16
#processors * * 2
read_data "cap.data"
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = lj"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

5
examples/electrode/check_hybrid/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0
#kspace_modify mesh 24 24 210
#kspace_modify order 5

36
examples/electrode/check_hybrid/ref.in Normal file
View File

@ -0,0 +1,36 @@
variable tag string ref
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
include "cgil.settings"
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = ref"
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

2
examples/electrode/check_intel/.gitignore vendored Normal file
View File

@ -0,0 +1,2 @@
log*
*csv

4766
examples/electrode/check_intel/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

23
examples/electrode/check_intel/cgil.settings Normal file
View File

@ -0,0 +1,23 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4
group carbon type 5
group carbon_neutral subtract carbon ele
neigh_modify exclude group carbon carbon_neutral

12
examples/electrode/check_intel/check.py Executable file
View File

@ -0,0 +1,12 @@
#!/bin/python3
import numpy as np
nointel = np.loadtxt("ref_nointel.csv", skiprows=1)
intel = np.loadtxt("ref_intel.csv", skiprows=1)
for matrix in [nointel, intel]:
assert matrix.shape == (6, 3), matrix.shape
diff0 = np.sum(np.abs(nointel - intel))
assert diff0 < 5e-4, diff0

17
examples/electrode/check_intel/check.sh Executable file
View File

@ -0,0 +1,17 @@
#!/bin/bash -e
lmpbin=$1
if [ ! -f $lmpbin ]; then
echo "LAMMPS binary '$lmpbin' is not a file"
exit 1
fi
rm -f *.csv
$lmpbin -i ref.in #> /dev/null
mv ref.csv ref_nointel.csv
$lmpbin -i ref.in -pk intel 0 omp 1 mode double -sf intel #> /dev/null
mv ref.csv ref_intel.csv
python3 check.py
#rm *.csv log.lammps*

5
examples/electrode/check_intel/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
kspace_modify gewald 0.21
kspace_modify slab 3.0
kspace_modify mesh 32 32 200
kspace_modify order 5

32
examples/electrode/check_intel/ref.in Normal file
View File

@ -0,0 +1,32 @@
variable tag string ref
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
include "p3m3dc.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
#processors * * 2
read_data "cap.data"
include "cgil.settings"
neigh_modify check no # force neighbor list rebuild for cost testing
#dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 symm on
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%20.15e) $(c_qbot:%20.15e)" file ${tag}.csv screen no title "#Fix electrode equal-style variable test setting = etypes"
thermo 10
timestep 2.5
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 50

26110
examples/electrode/example/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

30
examples/electrode/example/conp.in Normal file
View File

@ -0,0 +1,30 @@
atom_style full
units real
bond_style harmonic
angle_style harmonic
# Select one file for kspace settings
include "p3m3dc.settings"
# include "ew3dc.settings"
# include "ew2d.settings"
boundary p p f
pair_style lj/cut/coul/long 15 12
processors * * 2
read_data "cap.data"
include "mw.settings"
dump alldump all custom 1000 all.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 298.0 298.0 241
fix fxupdate bot electrode/conp -1.0 1.805132 couple top 1.0 write_inv fix_inv.csv symm on
compute q ele property/atom q
compute qtop top reduce sum c_q
compute qbot bot reduce sum c_q
compute qall ele reduce sum c_q
compute ctemp electrolyte temp
thermo 100
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 200

4
examples/electrode/example/ew2d.settings Normal file
View File

@ -0,0 +1,4 @@
kspace_style ewald/electrode 1.0e-6
kspace_modify gewald 0.2311815
kspace_modify kmax/ewald 10 10 29
kspace_modify slab ew2d

4
examples/electrode/example/ew3dc.settings Normal file
View File

@ -0,0 +1,4 @@
kspace_style ewald/electrode 1.0e-6
kspace_modify gewald 0.2311815
kspace_modify kmax/ewald 10 10 88
kspace_modify slab 3.0

17
examples/electrode/example/mw.settings Normal file
View File

@ -0,0 +1,17 @@
pair_modify mix arithmetic
pair_coeff 1 1 0.15540152963671128 3.166
pair_coeff 4 4 0.09999999999999999 2.584
pair_coeff 5 5 0.09999999999999999 4.401
pair_coeff 6 6 5.29 2.951
pair_coeff 7 7 5.29 2.951
pair_coeff 2 2 0.0 0.0
pair_coeff 3 3 0.0 0.0
bond_coeff 1 450 0.9999965664194801
bond_coeff 2 450 0.9999965664194801
angle_coeff 1 55 109.46937551578806
group SPC type 1 2 3
fix fx_rattle_SPC SPC shake 1.0e-4 20 0 b 1 2 a 1
group bot type 6
group top type 7
group ele type 6 7
group electrolyte type 1 2 3 4 5

5
examples/electrode/example/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-6
kspace_modify gewald 0.2311815
kspace_modify slab 3.0
kspace_modify mesh 24 24 210
kspace_modify order 5

4
examples/electrode/test_cgil_2d3d/.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
log.lammps*
fix_inv*
printout*
*lammpstrj

4766
examples/electrode/test_cgil_2d3d/cap.data Normal file
View File

File diff suppressed because it is too large Load Diff

19
examples/electrode/test_cgil_2d3d/cgil.settings Normal file
View File

@ -0,0 +1,19 @@
#Roy-Maroncelli BMI-PF6 + carbon
pair_modify mix arithmetic
pair_coeff 1 1 0.61185 4.38
pair_coeff 2 2 0.08604 3.41
pair_coeff 3 3 0.43738 5.04
pair_coeff 4 4 1.12572 5.06
pair_coeff 5 5 0.05497 3.37
bond_coeff 1 200 2.7076
bond_coeff 2 200 3.8213
angle_coeff 1 200 116.035
group BMI type 1 2 3
fix fx_rattle_BMI BMI shake 1e-4 20 0 b 1 2 a 1
group bot molecule 641
group top molecule 642
group ele union bot top
group electrolyte type 1 2 3 4

34
examples/electrode/test_cgil_2d3d/conp.in Normal file
View File

@ -0,0 +1,34 @@
variable tag getenv TAG
log log.lammps_${tag}
atom_style full
units real
bond_style harmonic
angle_style harmonic
# Select one file for kspace settings
if "${tag} == p3m" then "include p3m3dc.settings" &
elif "${tag} == ew3dc" "include ew3dc.settings" &
elif "${tag} == ew2d" "include ew2d.settings"
boundary p p f
pair_style lj/cut/coul/long 16 16
processors * * 2
read_data "cap.data"
include "cgil.settings"
dump alldump all custom 1000 all_${tag}.lammpstrj id mol type x y z q
fix fxnvt electrolyte nvt temp 500.0 500.0 100.
fix fxupdate bot electrode/conp -1.0 1.979 couple top 1.0 write_inv fix_inv_${tag}.csv symm on
#fix conp bot conp 1 top 1.979 2 log_conp pppm
variable q atom q
compute qtop top reduce sum v_q
compute qbot bot reduce sum v_q
compute qall ele reduce sum v_q
compute ctemp electrolyte temp
fix print all print 10 "$(step:%8.3g) $(c_qtop:%8.3e) $(c_qbot:%8.3e) $(c_ctemp:%8.3f) $(etotal:%12.3f)" file printout_${tag} screen no
thermo 10
thermo_style custom step c_ctemp epair etotal c_qtop c_qbot c_qall
run 500

4
examples/electrode/test_cgil_2d3d/ew2d.settings Normal file
View File

@ -0,0 +1,4 @@
kspace_style ewald/electrode 1.0e-6
kspace_modify gewald 0.2311815
kspace_modify kmax/ewald 10 10 29
kspace_modify slab ew2d

4
examples/electrode/test_cgil_2d3d/ew3dc.settings Normal file
View File

@ -0,0 +1,4 @@
kspace_style ewald/electrode 1.0e-6
#kspace_modify gewald 0.2311815
#kspace_modify kmax/ewald 10 10 88
kspace_modify slab 3.0

5
examples/electrode/test_cgil_2d3d/p3m3dc.settings Normal file
View File

@ -0,0 +1,5 @@
kspace_style pppm/electrode 1.0e-7
#kspace_modify gewald 0.2311815
kspace_modify slab 3.0
#kspace_modify mesh 24 24 210
#kspace_modify order 5

10
examples/electrode/test_cgil_2d3d/runtests.sh Normal file
View File

@ -0,0 +1,10 @@
LMP_EXEC=lmp_electrode
tags=(p3m ew3dc ew2d)
for t in ${tags[@]}
do
export TAG=$t
mpirun -np 4 $LMP_EXEC -i conp.in
done
paste printout_*

7
examples/electrode/thermo/.gitignore vendored Normal file
View File

@ -0,0 +1,7 @@
*.dat
*.xyz
log.lammps
mat.csv
final.inp
random.inp
nohup.out

131
examples/electrode/thermo/conp/control.inp Normal file
View File

@ -0,0 +1,131 @@
# Thermopotentiostat example calculation for nano-confined water
# between two electrodes (Ne) with applied voltage. Simulation
# can be performed with unmodified version of LAMMPS. Run as:
#
# lmp -in control.inp
#
# Output files:
# temp.dat temperature
# md.xyz trajectory
# final.inp final configuration for restarting
#
# F. Deissenbeck, C. Freysoldt, M. Todorova, J. Neugebauer, S. Wippermann,
# wippermann@mpie.de, 12.08.2020
units real
dimension 3
atom_style full
processors * * 2
kspace_style pppm/electrode 1.0e-5
kspace_modify slab 4.0 # non-periodic electrostatic boundary conditions in z-direction
kspace_modify gewald 0.2311815
boundary p p f
read_data ../sample.inp
include ../potential.inp
neighbor 2.0 bin
neigh_modify one 3000
# input parameters
# ----------------
variable dt equal 0.967553728 # time step [fs],
# 0.483776864 fs = 20 Hartree a.u.
variable dumpsteps equal 1000 # dump trajectory every ${dumpsteps} steps
group moving_particles type 1 2
group electrodes type 3 4
group Ne_left type 3
group Ne_right type 4
# NVE section (+ explicit thermostat if desired)
# ----------------------------------------------
timestep ${dt}
# init velocities to setup new claculation,
# not needed here since sample.inp is already equilibrated (at Phi0 = 0 V)
#velocity moving_particles create 350.0 4928459 rot yes dist gaussian
# apply explicit thermostat if desired,
# also not needed here since electrode charge distribution
# of thermopotentiostat has an explicit temperature (charge fluctuations)
# that actively controls the kinetic temperature
#fix 1 moving_particles temp/csvr 350.0 350.0 ${th_time} 54324
#fix 1 moving_particles langevin 350.0 350.0 ${th_time} 699483
# integrate Hamiltonian equations of motion -> pure NVE ensemble
fix 2 moving_particles nve
# ensure electrodes are not moving
fix freeze1 Ne_left setforce 0.0 0.0 0.0
fix freeze2 Ne_right setforce 0.0 0.0 0.0
# thermopotentiostat section
# --------------------------
fix conpid Ne_left electrode/conp 1.0 1.805132 couple Ne_right -1.0 symm on write_mat mat.csv
# output section
# --------------
dump 5 all custom ${dumpsteps} md.xyz id element xu yu zu q
dump_modify 5 format float %20.15g element H O Ne_a Ne_b sort id
# compute temperature
compute myTemp moving_particles temp
compute movppe moving_particles pe/atom
compute myPe moving_particles reduce sum c_movppe
compute myKe moving_particles ke
fix 6 moving_particles ave/time 1 1 1 c_myTemp c_myPe c_myKe file temp.dat format '%20.16f'
# compute charges
compute q electrodes property/atom q
compute q_left Ne_left reduce sum c_q
compute q_right Ne_right reduce sum c_q
compute q_all electrodes reduce sum c_q
thermo_style custom step etotal c_myTemp c_q_left c_q_right c_q_all f_conpid[1] f_conpid[2]
thermo 100
fix 3 moving_particles shake 0.0001 20 0 b 1 a 1
# run simulation and write final configuration for restart
# --------------------------------------------------------
# total number of timesteps
run 100000
# write final coordinates and velocities,
# nocoeff needed, or LAMMPS will not be able to restart from final.inp
write_data final.inp nocoeff

260
examples/electrode/thermo/flat/control.inp Normal file
View File

@ -0,0 +1,260 @@
# Thermopotentiostat example calculation for nano-confined water
# between two electrodes (Ne) with applied voltage. Simulation
# can be performed with unmodified version of LAMMPS. Run as:
#
# lmp -in control.inp
#
# Output files:
# temp.dat temperature
# traj_Phi.dat potential
# traj_n.dat charge
# md.xyz trajectory
# final.inp final configuration for restarting
#
# F. Deissenbeck, C. Freysoldt, M. Todorova, J. Neugebauer, S. Wippermann,
# wippermann@mpie.de, 12.08.2020
units real
dimension 3
atom_style full
kspace_style pppm/electrode 1.0e-5
kspace_modify slab 4.0 # non-periodic electrostatic boundary conditions in z-direction
kspace_modify gewald 0.2311815
boundary p p f
read_data ../sample.inp
include ../potential.inp
neighbor 2.0 bin
neigh_modify one 3000
# input parameters
# ----------------
variable dt equal 0.967553728 # time step [fs],
# 0.483776864 fs = 20 Hartree a.u.
variable th_temp equal 350.0 # temperature [K]
variable phi0 equal 2.0 # external applied potential [V]
variable po_time equal 100.0 # relaxation time potentiostat [fs]
variable dumpsteps equal 1000 # dump trajectory every ${dumpsteps} steps
# constants
# ---------
variable eps0 equal 0.005526349358057108 # [e-/(V Angstrom)]
variable kB equal 8.617333262145E-5 # [eV/K]
variable C equal v_eps0*v_Area/v_distance # [e-/V]
# setup - do not change (unless lowest electrode is not at z=2.0)
# ---------------------------------------------------------------
variable Area equal ($(xhi)-$(xlo))*($(yhi)-$(ylo)) # [A**2]
group moving_particles type 1 2
group electrodes type 3 4
group Ne_left type 3
group Ne_right type 4
variable Neatoms equal count(Ne_left)
# get electrode distance (to compute C) and
# one electrode atom (Ne) id (to update charge during simulation)
variable r atom gmask(Ne_right)*sqrt((z-10.0)*(z-10.0))
compute computedis all reduce max v_r # max dist between electrodes
compute ne_left_id Ne_left property/atom id
compute one_id_neleft Ne_left reduce min c_ne_left_id # get minimum id
thermo_style custom c_computedis c_one_id_neleft # dump to check
run 0
variable distance equal c_computedis
variable id_of_neleft equal c_one_id_neleft
# ----------------------- #
# | simulation protocol | #
# ----------------------- #
# compute dipole moments
# ----------------------
# total dipole
compute waterNe all chunk/atom bin/1d z lower $(zhi)
compute dip_tot all dipole/chunk waterNe
fix 5 all ave/time 1 1 1 c_dip_tot[*] file dipole_total.dat mode vector
# water-only dipole
compute onlywater moving_particles chunk/atom bin/1d z lower $(zhi)
compute dip_water all dipole/chunk onlywater
fix 43 all ave/time 1 1 1 c_dip_water[*] file dipole_water.dat mode vector
# electrode-only dipole, needed to compute electrode charge without recursion
compute elcharge electrodes chunk/atom bin/1d z lower $(zhi)
compute dip_elec all dipole/chunk elcharge
fix 83 all ave/time 1 1 1 c_dip_elec[*] file dipole_electrode.dat mode vector
# NVE section (+ explicit thermostat if desired)
# ----------------------------------------------
timestep ${dt}
# init velocities to setup new claculation,
# not needed here since sample.inp is already equilibrated (at Phi0 = 0 V)
#velocity moving_particles create 350.0 4928459 rot yes dist gaussian
# apply explicit thermostat if desired,
# also not needed here since electrode charge distribution
# of thermopotentiostat has an explicit temperature (charge fluctuations)
# that actively controls the kinetic temperature
# integrate Hamiltonian equations of motion -> pure NVE ensemble
fix 2 moving_particles nve
# ensure electrodes are not moving
fix freeze1 Ne_left setforce 0.0 0.0 0.0
fix freeze2 Ne_right setforce 0.0 0.0 0.0
# thermopotentiostat section
# --------------------------
# calculate electrode charge n from electrode-only dipole moment,
# done in this way in order to avoid a recursive loop due to lammps
# functional scripting language
variable charge_t1 equal -1*(c_dip_elec[1][3])/v_distance
# obtain voltage Phi from total dipole moment
# phi = -dip[3] / (eps0 * Area)
variable phi equal -1*c_dip_tot[1][3]/(v_eps0*v_Area)
# Ohm's Law potentiostat,
# dissipative term of Eq. 7 without fluctuation, cools system,
# only for testing purposes and by default deactivated
#variable phi_new equal v_phi0+(v_phi-v_phi0)*exp(-1*v_dt/v_po_time)
# CANONICAL SAMPLING BY STOCHASTIC THERMOPOTENTIOSTAT,
# equivalent to Eq. 7, just using Phi as the basic variable here,
# phi_new = phi0 + (phi - phi0) * np.exp(-dt/po_time) \
# + np.sqrt(kB*th_temp/C) * np.sqrt(1-np.exp(-2*dt/po_time))*np.random.normal(0.0,1.0),
# random number is read from file here to ensure reproducibility
# iterate random number every timestep with dummy fix
variable random file ../random.inp
variable phi_new equal v_phi0+(v_phi-v_phi0)*exp(-1*v_dt/v_po_time)+sqrt((v_kB*v_th_temp/v_C)*(1.0-exp(-2*v_dt/v_po_time)))*v_random
variable dummy equal next(random)
fix fxdummy all ave/time 1 1 1 v_dummy
# compute new charges on left and right electrodes
# dq = C * (phi_new - phi), q = q + dq
# performed specifically in this way to access the charge and hence phi_new only once,
# lammps scripting is functional, not imperative, and hence "reading" phi_new twice
# would also evaluate phi_new twice with different random numbers N and update charges on left electrode -> must avoid!
variable c equal (v_charge_t1+v_C*(v_phi_new-v_phi))
variable charge equal v_c/(1.0*v_Neatoms)
variable charge_2 equal -1.0*q[v_id_of_neleft]
# set charges on left and right electrodes,
# warning about unit cell being not charge-neutral can be ignored:
# after adapting left electrode first, unit cell is indeed non-neutral,
# after adapting right electrode subsequently it is neutral again
fix potentio_left Ne_left adapt 1 atom charge v_charge
fix potentio_right Ne_right adapt 1 atom charge v_charge_2
# output section
# --------------
dump 5 all custom ${dumpsteps} md.xyz id element xu yu zu q
dump_modify 5 format float %20.15g element H O Ne_a Ne_b sort id
fix potential all ave/time 1 1 1 v_phi file traj_Phi.dat
fix rho all ave/time 1 1 1 v_charge_t1 file traj_n.dat
# compute temperature
compute myTemp moving_particles temp
compute movppe moving_particles pe/atom
compute myPe moving_particles reduce sum c_movppe
compute myKe moving_particles ke
fix 6 moving_particles ave/time 1 1 1 c_myTemp c_myPe c_myKe file temp.dat format '%20.16f'
# compute charges
compute q electrodes property/atom q
compute q_left Ne_left reduce sum c_q
compute q_right Ne_right reduce sum c_q
compute q_all electrodes reduce sum c_q
thermo_style custom step etotal c_myTemp c_q_left c_q_right c_q_all
thermo 100
# compute temperature individually for H and O atoms if desired
#group movingH type 1
#group movingO type 2
#compute myTempH movingH temp
#compute myTempO movingO temp
#fix 7 moving_particles ave/time 1 1 1 c_myTempH c_myTempO file tempOH.dat format '%20.16f'
# rigid water => use shake-algorithm
# ----------------------------------
# must be at end of input file,
# if deactivating shake remember to choose smaller timestep,
# i.e. dt = 0.5 fs (or 20 Hartree atomic units) works fine
fix 3 moving_particles shake 0.0001 20 0 b 1 a 1
# run simulation and write final configuration for restart
# --------------------------------------------------------
# total number of timesteps
run 100000
# write final coordinates and velocities,
# nocoeff needed, or LAMMPS will not be able to restart from final.inp
write_data final.inp nocoeff

27
examples/electrode/thermo/potential.inp Normal file
View File

@ -0,0 +1,27 @@
# W.L. Jorgensen et.al., The Journal of Chemical Physics 79926 (1983);
# https://doi.org/10.1063/1.44586
# create groups
group H type 1
group O type 2
group Ne_a type 3
group Ne_b type 4
# set charges
set group O charge -0.830
set group H charge 0.415
set group Ne_a charge 0
set group Ne_b charge 0
# TIP3P Potential Parameters
pair_style lj/cut/coul/long 16.0
pair_coeff * * 0.0 0.0
pair_coeff 2 2 0.102 3.188
bond_style harmonic
bond_coeff 1 450 0.9572
angle_style harmonic
angle_coeff 1 55 104.52
# O - electrode (Ne) interaction
pair_coeff 2 3 0.102 3.188
pair_coeff 2 4 0.102 3.188

7
examples/electrode/thermo/rng.py Executable file
View File

@ -0,0 +1,7 @@
#!/usr/bin/env python3
import numpy as np
rng = np.random.default_rng(7)
vals = rng.standard_normal(10 ** 6 + 1)
np.savetxt("random.inp", vals)

14853
examples/electrode/thermo/sample.inp Normal file
View File

File diff suppressed because it is too large Load Diff

131
examples/electrode/thermo/thermo/control.inp Normal file
View File

@ -0,0 +1,131 @@
# Thermopotentiostat example calculation for nano-confined water
# between two electrodes (Ne) with applied voltage. Simulation
# can be performed with unmodified version of LAMMPS. Run as:
#
# lmp -in control.inp
#
# Output files:
# temp.dat temperature
# md.xyz trajectory
# final.inp final configuration for restarting
#
# F. Deissenbeck, C. Freysoldt, M. Todorova, J. Neugebauer, S. Wippermann,
# wippermann@mpie.de, 12.08.2020
units real
dimension 3
atom_style full
processors * * 2
kspace_style pppm/electrode 1.0e-5
kspace_modify slab 4.0 # non-periodic electrostatic boundary conditions in z-direction
kspace_modify gewald 0.2311815
boundary p p f
read_data ../sample.inp
include ../potential.inp
neighbor 2.0 bin
neigh_modify one 3000
# input parameters
# ----------------
variable dt equal 0.967553728 # time step [fs],
# 0.483776864 fs = 20 Hartree a.u.
variable dumpsteps equal 1000 # dump trajectory every ${dumpsteps} steps
group moving_particles type 1 2
group electrodes type 3 4
group Ne_left type 3
group Ne_right type 4
# NVE section (+ explicit thermostat if desired)
# ----------------------------------------------
timestep ${dt}
# init velocities to setup new claculation,
# not needed here since sample.inp is already equilibrated (at Phi0 = 0 V)
#velocity moving_particles create 350.0 4928459 rot yes dist gaussian
# apply explicit thermostat if desired,
# also not needed here since electrode charge distribution
# of thermopotentiostat has an explicit temperature (charge fluctuations)
# that actively controls the kinetic temperature
#fix 1 moving_particles temp/csvr 350.0 350.0 ${th_time} 54324
#fix 1 moving_particles langevin 350.0 350.0 ${th_time} 699483
# integrate Hamiltonian equations of motion -> pure NVE ensemble
fix 2 moving_particles nve
# ensure electrodes are not moving
fix freeze1 Ne_left setforce 0.0 0.0 0.0
fix freeze2 Ne_right setforce 0.0 0.0 0.0
# thermopotentiostat section
# --------------------------
fix thermoid Ne_left electrode/thermo 1.0 1.805132 couple Ne_right -1.0 temp 350 100 write_mat mat.csv
# output section
# --------------
dump 5 all custom ${dumpsteps} md.xyz id element xu yu zu q
dump_modify 5 format float %20.15g element H O Ne_a Ne_b sort id
# compute temperature
compute myTemp moving_particles temp
compute movppe moving_particles pe/atom
compute myPe moving_particles reduce sum c_movppe
compute myKe moving_particles ke
fix 6 moving_particles ave/time 1 1 1 c_myTemp c_myPe c_myKe file temp.dat format '%20.16f'
# compute charges
compute q electrodes property/atom q
compute q_left Ne_left reduce sum c_q
compute q_right Ne_right reduce sum c_q
compute q_all electrodes reduce sum c_q
thermo_style custom step etotal c_myTemp c_q_left c_q_right c_q_all f_thermoid[1] f_thermoid[2]
thermo 100
fix 3 moving_particles shake 0.0001 20 0 b 1 a 1
# run simulation and write final configuration for restart
# --------------------------------------------------------
# total number of timesteps
run 100000
# write final coordinates and velocities,
# nocoeff needed, or LAMMPS will not be able to restart from final.inp
write_data final.inp nocoeff

46
src/ELECTRODE/Install.sh Executable file
View File

@ -0,0 +1,46 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
# this is default Install.sh for all packages
# if package has an auxiliary library or a file with a dependency,
# then package dir has its own customized Install.sh
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*conp[^ \t]* //g' ../Makefile.package
sed -i -e 's|^PKG_LIB =[ \t]*|&-llinalg -L../../lib/linalg$(LIBOBJDIR) -lgfortran |' ../Makefile.package
fi
fi

95
src/ELECTRODE/boundary_correction.cpp Normal file
View File

@ -0,0 +1,95 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include "boundary_correction.h"
#include <iostream>
#include "atom.h"
#include "comm.h"
#include "force.h"
using namespace LAMMPS_NS;
using namespace std;
// use EW3DC slab correction
//
BoundaryCorrection::BoundaryCorrection(LAMMPS *lmp) : Pointers(lmp) {}
void BoundaryCorrection::setup(double x, double y, double z)
{
xprd_wire = x;
yprd_wire = y;
zprd_slab = z;
volume = x * y * z;
area = x * y;
qqrd2e = force->qqrd2e;
scale = 1.0;
}
void BoundaryCorrection::setup(double x, double y, double z, double g)
{
xprd_wire = x;
yprd_wire = y;
zprd_slab = z;
volume = x * y * z;
area = x * y;
qqrd2e = force->qqrd2e;
scale = 1.0;
g_ewald = g;
}
vector<int> BoundaryCorrection::gather_recvcounts(int n)
{
int const nprocs = comm->nprocs;
vector<int> recvcounts = vector<int>(nprocs);
MPI_Allgather(&n, 1, MPI_INT, &recvcounts.front(), 1, MPI_INT, world);
return recvcounts;
}
vector<int> BoundaryCorrection::gather_displs(vector<int> recvcounts)
{
int const nprocs = comm->nprocs;
vector<int> displs = vector<int>(nprocs);
displs[0] = 0;
for (int i = 1; i < nprocs; i++) displs[i] = displs[i - 1] + recvcounts[i - 1];
return displs;
}
vector<bigint> BoundaryCorrection::gather_jmat(bigint *imat)
{
int nlocal = atom->nlocal;
bigint ngroup = 0;
int ngrouplocal = 0;
for (int i = 0; i < nlocal; i++)
if (imat[i] > -1) ngrouplocal++;
MPI_Allreduce(&ngrouplocal, &ngroup, 1, MPI_INT, MPI_SUM, world);
vector<bigint> jmat_local = vector<bigint>(ngrouplocal);
for (int i = 0, n = 0; i < nlocal; i++) {
if (imat[i] < 0) continue;
jmat_local[n++] = imat[i];
}
// gather global matrix indexing
vector<bigint> jmat = vector<bigint>(ngroup);
vector<int> recvcounts = gather_recvcounts(ngrouplocal);
vector<int> displs = gather_displs(recvcounts);
MPI_Allgatherv(&jmat_local.front(), ngrouplocal, MPI_LMP_BIGINT, &jmat.front(),
&recvcounts.front(), &displs.front(), MPI_LMP_BIGINT, world);
return jmat;
}

49
src/ELECTRODE/boundary_correction.h Normal file
View File

@ -0,0 +1,49 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#ifndef LMP_BOUNDARY_CORRECTION_H
#define LMP_BOUNDARY_CORRECTION_H
#include "pointers.h"
namespace LAMMPS_NS {
class BoundaryCorrection : protected Pointers {
public:
BoundaryCorrection(LAMMPS *);
virtual void vector_corr(bigint *, double *){};
virtual void matrix_corr(bigint *, double **){};
virtual void compute_corr(double, int, int, double &, double *){};
void setup(double, double, double);
void setup(double, double, double, double);
protected:
double area;
double volume;
double xprd_wire;
double yprd_wire;
double zprd_slab;
double qqrd2e;
double scale;
double g_ewald;
std::vector<bigint> gather_jmat(bigint *);
std::vector<int> gather_recvcounts(int);
std::vector<int> gather_displs(std::vector<int>);
};
} // namespace LAMMPS_NS
#endif

20
src/ELECTRODE/electrode_accel_interface.cpp Normal file
View File

@ -0,0 +1,20 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include "electrode_accel_interface.h"
LAMMPS_NS::ElectrodeAccelInterface::~ElectrodeAccelInterface() {}

34
src/ELECTRODE/electrode_accel_interface.h Normal file
View File

@ -0,0 +1,34 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#ifndef LMP_ELECTRODE_ACCEL_INTERFACE_H
#define LMP_ELECTRODE_ACCEL_INTERFACE_H
#include "pointers.h"
namespace LAMMPS_NS {
class ElectrodeAccelInterface : protected Pointers {
public:
ElectrodeAccelInterface(class LAMMPS *lmp) : Pointers(lmp) {}
virtual ~ElectrodeAccelInterface();
virtual void intel_find_fix() {}
virtual void intel_pack_buffers() {}
};
} // namespace LAMMPS_NS
#endif

33
src/ELECTRODE/electrode_kspace.h Normal file
View File

@ -0,0 +1,33 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#ifndef LMP_ELECTRODE_KSPACE_H
#define LMP_ELECTRODE_KSPACE_H
#include "lmptype.h"
namespace LAMMPS_NS {
class ElectrodeKSpace {
public:
virtual void compute_vector(bigint *, double *) = 0;
virtual void compute_vector_corr(bigint *, double *) = 0;
virtual void compute_matrix(bigint *, double **) = 0;
virtual void compute_matrix_corr(bigint *, double **) = 0;
};
} // namespace LAMMPS_NS
#endif

201
src/ELECTRODE/electrode_matrix.cpp Normal file
View File

@ -0,0 +1,201 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include "electrode_matrix.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "kspace.h"
#include "math_const.h"
#include "neigh_list.h"
#include "pair.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
ElectrodeMatrix::ElectrodeMatrix(LAMMPS *lmp, int electrode_group, double eta) : Pointers(lmp)
{
igroup = electrode_group; // group of all electrode atoms
groupbit = group->bitmask[igroup];
ngroup = group->count(igroup);
this->eta = eta;
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::setup(std::map<tagint, int> tag_ids, class Pair *fix_pair,
class NeighList *fix_neighlist)
{
pair = fix_pair;
cutsq = pair->cutsq;
list = fix_neighlist;
electrode_kspace = dynamic_cast<ElectrodeKSpace *>(force->kspace);
if (electrode_kspace == nullptr) error->all(FLERR, "KSpace does not implement ElectrodeKSpace");
g_ewald = force->kspace->g_ewald;
tag_to_iele = tag_ids;
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::compute_array(double **array)
{
// setting all entries of coulomb matrix to zero
size_t nbytes = sizeof(double) * ngroup * ngroup;
if (nbytes) memset(&array[0][0], 0, nbytes);
MPI_Barrier(world);
double kspace_time = MPI_Wtime();
update_mpos();
electrode_kspace->compute_matrix(&mpos[0], array);
MPI_Barrier(world);
if (comm->me == 0)
utils::logmesg(lmp, fmt::format("KSpace time: {}\n", MPI_Wtime() - kspace_time));
pair_contribution(array);
self_contribution(array);
electrode_kspace->compute_matrix_corr(&mpos[0], array);
// reduce coulomb matrix with contributions from all procs
// all procs need to know full matrix for matrix inversion
for (int i = 0; i < ngroup; i++) {
MPI_Allreduce(MPI_IN_PLACE, &array[i][0], ngroup, MPI_DOUBLE, MPI_SUM, world);
}
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::pair_contribution(double **array)
{
int inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz;
double r, rinv, rsq, grij, etarij, expm2, t, erfc, aij;
int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x;
tagint *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double etaij = eta * eta / sqrt(2.0 * eta * eta); // see mw ewald theory eq. (29)-(30)
// neighbor list will be ready because called from post_neighbor
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
// skip if I,J are not in 2 groups
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
// skip if atom I is not in either group
if (!(mask[i] & groupbit)) continue;
bigint const ipos = mpos[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
// real-space part of matrix is symmetric
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
delx = xtmp - x[j][0]; // neighlists take care of pbc
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
rinv = 1.0 / r;
aij = rinv;
grij = g_ewald * r;
expm2 = exp(-grij * grij);
t = 1.0 / (1.0 + EWALD_P * grij);
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
aij *= erfc;
etarij = etaij * r;
expm2 = exp(-etarij * etarij);
t = 1.0 / (1.0 + EWALD_P * etarij);
erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
aij -= erfc * rinv;
// newton on or off?
if (!(newton_pair || j < nlocal)) aij *= 0.5;
bigint jpos = tag_to_iele[tag[j]];
array[ipos][jpos] += aij;
array[jpos][ipos] += aij;
}
}
}
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::self_contribution(double **array)
{
int nlocal = atom->nlocal;
int *mask = atom->mask;
const double selfint = 2.0 / MY_PIS * g_ewald;
const double preta = MY_SQRT2 / MY_PIS;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { array[mpos[i]][mpos[i]] += preta * eta - selfint; }
}
/* ---------------------------------------------------------------------- */
void ElectrodeMatrix::update_mpos()
{
int const nall = atom->nlocal + atom->nghost;
int *tag = atom->tag;
int *mask = atom->mask;
mpos = std::vector<bigint>(nall, -1);
for (int i = 0; i < nall; i++) {
if (mask[i] & groupbit)
mpos[i] = tag_to_iele[tag[i]];
else
mpos[i] = -1;
}
}

91
src/ELECTRODE/electrode_matrix.h Normal file
View File

@ -0,0 +1,91 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include <map>
#include "electrode_kspace.h"
#include "pointers.h"
namespace LAMMPS_NS {
class ElectrodeMatrix : protected Pointers {
public:
ElectrodeMatrix(class LAMMPS *, int, double);
~ElectrodeMatrix() {}
void setup(std::map<tagint, int>, class Pair *, class NeighList *);
void compute_array(double **);
int igroup;
private:
int groupbit;
bigint ngroup;
double **cutsq;
double g_ewald, eta;
std::map<tagint, int> tag_to_iele;
bool assigned;
std::vector<bigint> mpos;
class Pair *pair;
class NeighList *list;
class ElectrodeKSpace *electrode_kspace;
void update_mpos();
void pair_contribution(double **);
void self_contribution(double **);
double calc_erfc(double);
};
} // namespace LAMMPS_NS
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute group/group group ID does not exist
Self-explanatory.
E: Compute group/group molecule requires molecule IDs
UNDOCUMENTED
E: No pair style defined for compute group/group
Cannot calculate group interactions without a pair style defined.
E: Pair style does not support compute group/group
The pair_style does not have a single() function, so it cannot be
invoked by the compute group/group command.
E: No KSpace style defined for compute group/group
Self-explanatory.
E: KSpace style does not support compute group/group
Self-explanatory.
W: Both groups in compute group/group have a net charge; the KSpace boundary
correction to energy will be non-zero
Self-explanatory.
*/

206
src/ELECTRODE/electrode_vector.cpp Normal file
View File

@ -0,0 +1,206 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include "electrode_vector.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "kspace.h"
#include "neigh_list.h"
#include "pair.h"
using namespace LAMMPS_NS;
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
ElectrodeVector::ElectrodeVector(LAMMPS *lmp, int electrode_group, double eta) : Pointers(lmp)
{
igroup = electrode_group; // group of all electrode atoms
groupbit = group->bitmask[igroup];
ngroup = group->count(igroup);
vector = new double[ngroup](); // init to zero
this->eta = eta;
setup_time_total = 0;
reduce_time_total = 0;
kspace_time_total = 0;
pair_time_total = 0;
boundary_time_total = 0;
b_time_total = 0;
alloc_time_total = 0;
mpos_time_total = 0;
}
/* ---------------------------------------------------------------------- */
ElectrodeVector::~ElectrodeVector()
{
if (comm->me == 0) {
utils::logmesg(lmp, fmt::format("B time: {}\n", b_time_total));
utils::logmesg(lmp, fmt::format("B kspace time: {}\n", kspace_time_total));
utils::logmesg(lmp, fmt::format("B pair time: {}\n", pair_time_total));
utils::logmesg(lmp, fmt::format("B boundary time: {}\n", boundary_time_total));
utils::logmesg(lmp, fmt::format("B setup time: {}\n", setup_time_total));
utils::logmesg(lmp, fmt::format("B reduce time: {}\n", reduce_time_total));
utils::logmesg(lmp, fmt::format("B alloc time: {}\n", alloc_time_total));
utils::logmesg(lmp, fmt::format("B mpos time: {}\n", mpos_time_total));
}
delete[] vector;
}
/* ---------------------------------------------------------------------- */
void ElectrodeVector::setup(std::map<tagint, int> tag_ids, class Pair *fix_pair,
class NeighList *fix_neighlist)
{
pair = fix_pair;
cutsq = pair->cutsq;
list = fix_neighlist;
electrode_kspace = dynamic_cast<ElectrodeKSpace *>(force->kspace);
if (electrode_kspace == nullptr) error->all(FLERR, "KSpace does not implement ElectrodeKSpace");
g_ewald = force->kspace->g_ewald;
tag_to_iele = tag_ids;
}
/* ---------------------------------------------------------------------- */
void ElectrodeVector::compute_vector()
{
MPI_Barrier(world);
double start_time = MPI_Wtime();
// setup
double setup_start_time = MPI_Wtime();
update_mpos();
for (int i = 0; i < ngroup; i++) vector[i] = 0.;
MPI_Barrier(world);
setup_time_total += MPI_Wtime() - setup_start_time;
// pair
double pair_start_time = MPI_Wtime();
pair_contribution();
MPI_Barrier(world);
pair_time_total += MPI_Wtime() - pair_start_time;
// kspace
double kspace_start_time = MPI_Wtime();
electrode_kspace->compute_vector(&mpos[0], vector);
MPI_Barrier(world);
kspace_time_total += MPI_Wtime() - kspace_start_time;
// boundary
double boundary_start_time = MPI_Wtime();
electrode_kspace->compute_vector_corr(&mpos[0], vector);
MPI_Barrier(world);
boundary_time_total += MPI_Wtime() - boundary_start_time;
// reduce
double reduce_start_time = MPI_Wtime();
MPI_Allreduce(MPI_IN_PLACE, vector, ngroup, MPI_DOUBLE, MPI_SUM, world);
MPI_Barrier(world);
reduce_time_total += MPI_Wtime() - reduce_start_time;
b_time_total += MPI_Wtime() - start_time;
}
/* ---------------------------------------------------------------------- */
void ElectrodeVector::pair_contribution()
{
double **x = atom->x;
double *q = atom->q;
int *type = atom->type;
int *mask = atom->mask;
// neighbor list will be ready because called from post_neighbor
int const nlocal = atom->nlocal;
int const inum = list->inum;
int *ilist = list->ilist;
int *numneigh = list->numneigh;
int **firstneigh = list->firstneigh;
int newton_pair = force->newton_pair;
for (int ii = 0; ii < inum; ii++) {
int const i = ilist[ii];
bool const i_in_electrode = (mask[i] & groupbit);
double const xtmp = x[i][0];
double const ytmp = x[i][1];
double const ztmp = x[i][2];
int itype = type[i];
int *jlist = firstneigh[i];
int jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
int const j = jlist[jj] & NEIGHMASK;
bool const j_in_electrode = (mask[j] & groupbit);
if (i_in_electrode == j_in_electrode) continue;
double const delx = xtmp - x[j][0]; // neighlists take care of pbc
double const dely = ytmp - x[j][1];
double const delz = ztmp - x[j][2];
double const rsq = delx * delx + dely * dely + delz * delz;
int jtype = type[j];
if (rsq >= cutsq[itype][jtype]) continue;
double const r = sqrt(rsq);
double const rinv = 1.0 / r;
double aij = rinv;
aij *= calc_erfc(g_ewald * r);
aij -= calc_erfc(eta * r) * rinv;
if (!(newton_pair || j < nlocal)) aij *= 0.5;
if (i_in_electrode) {
vector[mpos[i]] += aij * q[j];
} else if (j_in_electrode) {
vector[mpos[j]] += aij * q[i];
}
}
}
}
/* ---------------------------------------------------------------------- */
void ElectrodeVector::update_mpos()
{
MPI_Barrier(world);
double alloc_start = MPI_Wtime();
int const nall = atom->nlocal + atom->nghost;
int *tag = atom->tag;
int *mask = atom->mask;
mpos = std::vector<bigint>(nall, -1);
MPI_Barrier(world);
alloc_time_total += MPI_Wtime() - alloc_start;
double mpos_start = MPI_Wtime();
for (int i = 0; i < nall; i++) {
if (mask[i] & groupbit)
mpos[i] = tag_to_iele[tag[i]];
else
mpos[i] = -1;
}
MPI_Barrier(world);
mpos_time_total += MPI_Wtime() - mpos_start;
}
/* ---------------------------------------------------------------------- */
double ElectrodeVector::calc_erfc(double x)
{
double expm2 = exp(-x * x);
double t = 1.0 / (1.0 + EWALD_P * x);
return t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
}

61
src/ELECTRODE/electrode_vector.h Normal file
View File

@ -0,0 +1,61 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://lammps.sandia.gov/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Ludwig Ahrens-Iwers (TUHH), Shern Tee (UQ), Robert Meißner (TUHH)
------------------------------------------------------------------------- */
#include <map>
#include "electrode_kspace.h"
#include "pointers.h"
namespace LAMMPS_NS {
class ElectrodeVector : protected Pointers {
public:
ElectrodeVector(class LAMMPS *, int, double);
~ElectrodeVector();
void setup(std::map<tagint, int>, class Pair *, class NeighList *);
void compute_vector();
double *vector;
int igroup;
private:
int groupbit;
bigint ngroup;
double **cutsq;
double g_ewald, eta;
std::map<tagint, int> tag_to_iele;
std::vector<bigint> mpos;
class Pair *pair;
class NeighList *list;
class ElectrodeKSpace *electrode_kspace;
void update_mpos();
void pair_contribution();
double calc_erfc(double);
double setup_time_total;
double reduce_time_total;
double kspace_time_total;
double pair_time_total;
double boundary_time_total;
double b_time_total;
double alloc_time_total;
double mpos_time_total;
};
} // namespace LAMMPS_NS

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