implement compute_array

This commit is contained in:
srtee
2022-05-02 23:24:00 +10:00
parent eb5d867adf
commit 0855d62bfd
9 changed files with 126 additions and 259 deletions

View File

@ -130,48 +130,49 @@ file, for code developers to track optimization.
fix_modify ID timer on/off
The *fix_modify set* options allow calculated quantities to be accessed via
internal variables. Currently four types of quantities can be accessed:
----------
.. code-block:: LAMMPS
These fixes compute a global array which can be accessed by various
:doc:`output commands <Howto_output>`. The array has *N* rows and *2N+2*
columns, where *N* is the number of electrode groups managed by the fix. For the
*I*-th row of the array, the elements are:
fix-modify ID set v group-ID variablename
fix-modify ID set qsb group-ID variablename
fix-modify ID set mc group-ID1 group-ID2 variablename
fix-modify ID set me group-ID1 group-ID2 variablename
* array[I][1] = potential instantaneously applied to group *I*
* array[I][1] = total charge that group *I* would have had *if it were at 0 V
applied potential*
* array[I][2 to *N* + 1] = the *N* entries of the *I*-th row of the electrode
capacitance matrix (definition follows)
* array[I][*N* + 2 to *2N* + 1] = the *N* entries of the *I*-th row of the electrode
elastance matrix (the inverse of the electrode capacitance matrix)
One use case is to output the potential that is internally calculated and
applied to each electrode group by *fix electrode/conq* or *fix electrode/thermo*.
For that case the *v* option makes *fix electrode* update the variable
*variablename* with the potential applied to group *group-ID*, where *group-ID*
must be a group whose charges are updated by *fix electrode* and *variablename*
must be an internal-style variable:
The first column is handy for *fix electrode/conq* and *fix
electrode/thermo* simulations, to output the potentials dynamically calculated and
imposed by the fix (see following formula). The potential is in units of electric field times
distance, which gives volts for most (but not all) :doc:`units <units>` settings.
.. code-block:: LAMMPS
fix conq bot electrode/conq -1.0 1.979 couple top 1.0
variable vbot internal 0.0
fix_modify conq set v bot vbot
The *qsb* option similarly outputs the total updated charge of the group if its
potential were 0.0V. The *mc* option requires two *group-IDs*, and outputs the
entry \{*group-ID1*, *group-ID2*\} of the (symmetric) *macro-capacitance* matrix
(MC) which relates the electrodes' applied potentials (V), total charges (Q), and
total charges at 0.0 V (Qsb):
The :math:`N \times N` electrode capacitance matrix, denoted :math:`\mathbf{C}`
in the following equation, summarizes how the total charge induced on each electrode
(:math:`\mathbf{Q}` as an *N*-vector) is related to the potential on each
electrode, :math:`\mathbf{V}`, and the charge-at-0V :math:`\mathbf{Q}_{0V}`
(which is influenced by the local electrolyte structure):
.. math::
\mathbf{Q} = \mathbf{Q}_{SB} + \mathbf{MC} \cdot \mathbf{V}
\mathbf{Q} = \mathbf{Q}_{0V} + \mathbf{C} \cdot \mathbf{V}
Lastly, the *me* option also requires two *group-IDs* and outputs the entry
\{*group-ID1*, *group-ID2*\} of the *macro-elastance* matrix, which is the
inverse of the macro-capacitance matrix. (As the names denote, the
macro-capacitance matrix gives electrode charges from potentials, and the
macro-elastance matrix gives electrode potentials from charges).
The charge-at-0V, electrode capacitance and elastance matrices are internally used to
calculate the potentials required to induce the specified total electrode charges
in *fix electrode/conq* and *fix electrode/thermo*. With the *symm on* option,
the electrode capacitance matrix would be singular, and thus its last row is
replaced with *N* copies of its top-left entry (:math:`\mathbf{C}_{11}`) for
invertibility.
.. warning::
----------
Positions of electrode particles have to be immobilized at all times.
Restrictions
""""""""""""
Positions of electrode particles have to be immobilized at all times.
The parallelization for the fix works best if electrode atoms are evenly
distributed across processors. For a system with two electrodes at the bottom

View File

@ -7,12 +7,6 @@ kspace_modify slab 3.0
fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes 5 # conq doesn't take symm option
# ask fix conq to output electrode potentials to internal variables
variable vbot internal 0.0
variable vtop internal 0.0
fix_modify conq set v bot vbot
fix_modify conq set v top vtop
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1][1] f_conq[2][1]
run 500

View File

@ -7,28 +7,14 @@ kspace_modify slab 3.0
fix conp bot electrode/conp v_vbot 1.979 couple top v_vtop etypes 5
# get the four entries of electrode elastance matrix
variable me00 internal 0.0
variable me01 internal 0.0
variable me10 internal 0.0
variable me11 internal 0.0
fix_modify conp set me bot bot me00
fix_modify conp set me bot top me01
fix_modify conp set me top bot me10
fix_modify conp set me top top me11
# get the 0V charges (qsb), and excess charge required to reach preset total charges
variable qsb_bot internal 0.0
variable qsb_top internal 0.0
fix_modify conp set qsb bot qsb_bot
fix_modify conp set qsb top qsb_top
variable qex_bot equal -1.0-v_qsb_bot # difference between desired and 0V charge
variable qex_top equal 1.0-v_qsb_top # difference between desired and 0V charge
variable qex_bot equal -1.0-f_conp[1][2] # difference between desired and 0V charge
variable qex_top equal 1.0-f_conp[2][2] # difference between desired and 0V charge
# calculate imposed potential as elastance * excess charge
# note: fix will wait until the run setup to look for its potential variables
variable vbot equal v_me00*v_qex_bot+v_me01*v_qex_top
variable vtop equal v_me10*v_qex_bot+v_me11*v_qex_top
# which is why we can define variable names *after* fix conp without error
variable vbot equal f_conp[1][5]*v_qex_bot+f_conp[1][6]*v_qex_top
variable vtop equal f_conp[2][5]*v_qex_bot+f_conp[2][6]*v_qex_top
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop

View File

@ -67,14 +67,8 @@ kspace_modify slab 3.0
fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes 5 # conq doesn't take symm option
832 atoms in group conp_group
# ask fix conq to output electrode potentials to internal variables
variable vbot internal 0.0
variable vtop internal 0.0
fix_modify conq set v bot vbot
fix_modify conq set v top vtop
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1][1] f_conq[2][1]
run 500
PPPM/electrode initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
@ -109,7 +103,7 @@ Neighbor list info ...
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 56.89 | 56.89 | 56.89 Mbytes
Step Temp c_ctemp E_pair TotEng c_qbot c_qtop v_vbot v_vtop
Step Temp c_ctemp E_pair TotEng c_qbot c_qtop f_conq[1][1] f_conq[2][1]
0 0 0 25136984 25136984 -1 1 -9.931852 10.097344
50 20.206425 72.797911 25136825 25137033 -1 1 -9.4359366 9.5964514
100 55.931663 201.50563 25136587 25137163 -1 1 -8.0440112 8.1861787
@ -121,22 +115,22 @@ Per MPI rank memory allocation (min/avg/max) = 56.89 | 56.89 | 56.89 Mbytes
400 122.8443 442.57252 25136991 25138256 -1 1 -3.4311003 3.3941657
450 128.63713 463.44243 25137048 25138373 -1 1 -4.132871 3.9852959
500 131.18361 472.61665 25137142 25138493 -1 1 -4.5104095 4.2567261
Loop time of 48.9361 on 1 procs for 500 steps with 3776 atoms
Loop time of 42.3981 on 1 procs for 500 steps with 3776 atoms
Performance: 0.883 ns/day, 27.187 hours/ns, 10.217 timesteps/s
393.9% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 1.019 ns/day, 23.555 hours/ns, 11.793 timesteps/s
393.2% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 17.652 | 17.652 | 17.652 | 0.0 | 36.07
Bond | 0.0010418 | 0.0010418 | 0.0010418 | 0.0 | 0.00
Kspace | 16.566 | 16.566 | 16.566 | 0.0 | 33.85
Neigh | 0.21584 | 0.21584 | 0.21584 | 0.0 | 0.44
Comm | 0.04167 | 0.04167 | 0.04167 | 0.0 | 0.09
Output | 0.0014585 | 0.0014585 | 0.0014585 | 0.0 | 0.00
Modify | 14.445 | 14.445 | 14.445 | 0.0 | 29.52
Other | | 0.0134 | | | 0.03
Pair | 15.337 | 15.337 | 15.337 | 0.0 | 36.17
Bond | 0.00098464 | 0.00098464 | 0.00098464 | 0.0 | 0.00
Kspace | 16.11 | 16.11 | 16.11 | 0.0 | 38.00
Neigh | 0.19136 | 0.19136 | 0.19136 | 0.0 | 0.45
Comm | 0.039778 | 0.039778 | 0.039778 | 0.0 | 0.09
Output | 0.0015193 | 0.0015193 | 0.0015193 | 0.0 | 0.00
Modify | 10.705 | 10.705 | 10.705 | 0.0 | 25.25
Other | | 0.01253 | | | 0.03
Nlocal: 3776 ave 3776 max 3776 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -150,4 +144,4 @@ Ave neighs/atom = 456.98835
Ave special neighs/atom = 0.50847458
Neighbor list builds = 6
Dangerous builds = 0
Total wall time: 0:01:43
Total wall time: 0:01:34

View File

@ -39,7 +39,7 @@ Finding 1-2 1-3 1-4 neighbors ...
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.020 seconds
read_data CPU = 0.018 seconds
group bot molecule 641
416 atoms in group bot
@ -68,14 +68,8 @@ kspace_modify slab 3.0
fix conq bot electrode/conq -1.0 1.979 couple top 1.0 etypes 5 # conq doesn't take symm option
832 atoms in group conp_group
# ask fix conq to output electrode potentials to internal variables
variable vbot internal 0.0
variable vtop internal 0.0
fix_modify conq set v bot vbot
fix_modify conq set v top vtop
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop f_conq[1][1] f_conq[2][1]
run 500
PPPM/electrode initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
@ -110,7 +104,7 @@ Neighbor list info ...
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 23.63 | 27.46 | 31.29 Mbytes
Step Temp c_ctemp E_pair TotEng c_qbot c_qtop v_vbot v_vtop
Step Temp c_ctemp E_pair TotEng c_qbot c_qtop f_conq[1][1] f_conq[2][1]
0 0 0 25136984 25136984 -1 1 -9.931852 10.097344
50 20.206425 72.797911 25136825 25137033 -1 1 -9.4359366 9.5964514
100 55.931663 201.50563 25136587 25137163 -1 1 -8.0440112 8.1861787
@ -122,22 +116,22 @@ Per MPI rank memory allocation (min/avg/max) = 23.63 | 27.46 | 31.29 Mbytes
400 122.8443 442.57252 25136991 25138256 -1 1 -3.4311003 3.3941657
450 128.63713 463.44243 25137048 25138373 -1 1 -4.132871 3.9852959
500 131.18361 472.61665 25137142 25138493 -1 1 -4.5104095 4.2567261
Loop time of 28.8336 on 4 procs for 500 steps with 3776 atoms
Loop time of 23.02 on 4 procs for 500 steps with 3776 atoms
Performance: 1.498 ns/day, 16.019 hours/ns, 17.341 timesteps/s
94.1% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 1.877 ns/day, 12.789 hours/ns, 21.720 timesteps/s
93.3% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.7721 | 5.9353 | 6.144 | 6.0 | 20.58
Bond | 0.00057855 | 0.00067043 | 0.00074793 | 0.0 | 0.00
Kspace | 13.485 | 13.694 | 13.857 | 4.0 | 47.49
Neigh | 0.092021 | 0.092044 | 0.092068 | 0.0 | 0.32
Comm | 0.11486 | 0.11638 | 0.11801 | 0.4 | 0.40
Output | 0.00090452 | 0.001109 | 0.0017097 | 1.0 | 0.00
Modify | 8.974 | 8.9761 | 8.978 | 0.1 | 31.13
Other | | 0.01837 | | | 0.06
Pair | 4.0178 | 4.1947 | 4.3354 | 6.8 | 18.22
Bond | 0.00061117 | 0.00068637 | 0.00072249 | 0.0 | 0.00
Kspace | 11.853 | 11.994 | 12.171 | 4.0 | 52.10
Neigh | 0.060425 | 0.060489 | 0.06052 | 0.0 | 0.26
Comm | 0.091334 | 0.092251 | 0.093191 | 0.3 | 0.40
Output | 0.00067425 | 0.00089541 | 0.0014804 | 0.0 | 0.00
Modify | 6.6612 | 6.662 | 6.6627 | 0.0 | 28.94
Other | | 0.01505 | | | 0.07
Nlocal: 944 ave 948 max 940 min
Histogram: 1 0 0 1 0 0 1 0 0 1
@ -151,4 +145,4 @@ Ave neighs/atom = 456.98835
Ave special neighs/atom = 0.50847458
Neighbor list builds = 6
Dangerous builds = 0
Total wall time: 0:00:51
Total wall time: 0:00:44

View File

@ -38,7 +38,7 @@ Finding 1-2 1-3 1-4 neighbors ...
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.031 seconds
read_data CPU = 0.012 seconds
group bot molecule 641
416 atoms in group bot
@ -67,28 +67,14 @@ kspace_modify slab 3.0
fix conp bot electrode/conp v_vbot 1.979 couple top v_vtop etypes 5
832 atoms in group conp_group
# get the four entries of electrode elastance matrix
variable me00 internal 0.0
variable me01 internal 0.0
variable me10 internal 0.0
variable me11 internal 0.0
fix_modify conp set me bot bot me00
fix_modify conp set me bot top me01
fix_modify conp set me top bot me10
fix_modify conp set me top top me11
# get the 0V charges (qsb), and excess charge required to reach preset total charges
variable qsb_bot internal 0.0
variable qsb_top internal 0.0
fix_modify conp set qsb bot qsb_bot
fix_modify conp set qsb top qsb_top
variable qex_bot equal -1.0-v_qsb_bot # difference between desired and 0V charge
variable qex_top equal 1.0-v_qsb_top # difference between desired and 0V charge
variable qex_bot equal -1.0-f_conp[1][2] # difference between desired and 0V charge
variable qex_top equal 1.0-f_conp[2][2] # difference between desired and 0V charge
# calculate imposed potential as elastance * excess charge
# note: fix will wait until the run setup to look for its potential variables
variable vbot equal v_me00*v_qex_bot+v_me01*v_qex_top
variable vtop equal v_me10*v_qex_bot+v_me11*v_qex_top
# which is why we can define variable names *after* fix conp without error
variable vbot equal f_conp[1][5]*v_qex_bot+f_conp[1][6]*v_qex_top
variable vtop equal f_conp[2][5]*v_qex_bot+f_conp[2][6]*v_qex_top
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
@ -138,22 +124,22 @@ Per MPI rank memory allocation (min/avg/max) = 56.89 | 56.89 | 56.89 Mbytes
400 122.8443 442.57252 25136991 25138256 -1 1 -3.4311003 3.3941657
450 128.63713 463.44243 25137048 25138373 -1 1 -4.132871 3.9852959
500 131.18361 472.61665 25137142 25138493 -1 1 -4.5104095 4.2567261
Loop time of 62.9692 on 1 procs for 500 steps with 3776 atoms
Loop time of 42.7942 on 1 procs for 500 steps with 3776 atoms
Performance: 0.686 ns/day, 34.983 hours/ns, 7.940 timesteps/s
393.7% CPU use with 1 MPI tasks x no OpenMP threads
Performance: 1.009 ns/day, 23.775 hours/ns, 11.684 timesteps/s
394.4% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 20.307 | 20.307 | 20.307 | 0.0 | 32.25
Bond | 0.0020074 | 0.0020074 | 0.0020074 | 0.0 | 0.00
Kspace | 23.562 | 23.562 | 23.562 | 0.0 | 37.42
Neigh | 0.26149 | 0.26149 | 0.26149 | 0.0 | 0.42
Comm | 0.059436 | 0.059436 | 0.059436 | 0.0 | 0.09
Output | 0.0023888 | 0.0023888 | 0.0023888 | 0.0 | 0.00
Modify | 18.756 | 18.756 | 18.756 | 0.0 | 29.79
Other | | 0.01897 | | | 0.03
Pair | 15.743 | 15.743 | 15.743 | 0.0 | 36.79
Bond | 0.00097932 | 0.00097932 | 0.00097932 | 0.0 | 0.00
Kspace | 16.028 | 16.028 | 16.028 | 0.0 | 37.45
Neigh | 0.21147 | 0.21147 | 0.21147 | 0.0 | 0.49
Comm | 0.040091 | 0.040091 | 0.040091 | 0.0 | 0.09
Output | 0.001878 | 0.001878 | 0.001878 | 0.0 | 0.00
Modify | 10.756 | 10.756 | 10.756 | 0.0 | 25.13
Other | | 0.01243 | | | 0.03
Nlocal: 3776 ave 3776 max 3776 min
Histogram: 1 0 0 0 0 0 0 0 0 0
@ -167,4 +153,4 @@ Ave neighs/atom = 456.98835
Ave special neighs/atom = 0.50847458
Neighbor list builds = 6
Dangerous builds = 0
Total wall time: 0:02:26
Total wall time: 0:01:35

View File

@ -68,28 +68,14 @@ kspace_modify slab 3.0
fix conp bot electrode/conp v_vbot 1.979 couple top v_vtop etypes 5
832 atoms in group conp_group
# get the four entries of electrode elastance matrix
variable me00 internal 0.0
variable me01 internal 0.0
variable me10 internal 0.0
variable me11 internal 0.0
fix_modify conp set me bot bot me00
fix_modify conp set me bot top me01
fix_modify conp set me top bot me10
fix_modify conp set me top top me11
# get the 0V charges (qsb), and excess charge required to reach preset total charges
variable qsb_bot internal 0.0
variable qsb_top internal 0.0
fix_modify conp set qsb bot qsb_bot
fix_modify conp set qsb top qsb_top
variable qex_bot equal -1.0-v_qsb_bot # difference between desired and 0V charge
variable qex_top equal 1.0-v_qsb_top # difference between desired and 0V charge
variable qex_bot equal -1.0-f_conp[1][2] # difference between desired and 0V charge
variable qex_top equal 1.0-f_conp[2][2] # difference between desired and 0V charge
# calculate imposed potential as elastance * excess charge
# note: fix will wait until the run setup to look for its potential variables
variable vbot equal v_me00*v_qex_bot+v_me01*v_qex_top
variable vtop equal v_me10*v_qex_bot+v_me11*v_qex_top
# which is why we can define variable names *after* fix conp without error
variable vbot equal f_conp[1][5]*v_qex_bot+f_conp[1][6]*v_qex_top
variable vtop equal f_conp[2][5]*v_qex_bot+f_conp[2][6]*v_qex_top
thermo 50
thermo_style custom step temp c_ctemp epair etotal c_qbot c_qtop v_vbot v_vtop
@ -139,22 +125,22 @@ Per MPI rank memory allocation (min/avg/max) = 23.63 | 27.46 | 31.29 Mbytes
400 122.8443 442.57252 25136991 25138256 -1 1 -3.4311003 3.3941657
450 128.63713 463.44243 25137048 25138373 -1 1 -4.132871 3.9852959
500 131.18361 472.61665 25137142 25138493 -1 1 -4.5104095 4.2567261
Loop time of 33.4031 on 4 procs for 500 steps with 3776 atoms
Loop time of 23.4395 on 4 procs for 500 steps with 3776 atoms
Performance: 1.293 ns/day, 18.557 hours/ns, 14.969 timesteps/s
94.3% CPU use with 4 MPI tasks x no OpenMP threads
Performance: 1.843 ns/day, 13.022 hours/ns, 21.331 timesteps/s
93.3% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 7.1262 | 7.3913 | 7.611 | 6.8 | 22.13
Bond | 0.0007191 | 0.00079089 | 0.00087005 | 0.0 | 0.00
Kspace | 15.139 | 15.358 | 15.623 | 4.7 | 45.98
Neigh | 0.10374 | 0.10377 | 0.10383 | 0.0 | 0.31
Comm | 0.14245 | 0.14353 | 0.14563 | 0.3 | 0.43
Output | 0.0012987 | 0.0015671 | 0.0022434 | 1.0 | 0.00
Modify | 10.381 | 10.383 | 10.384 | 0.0 | 31.08
Other | | 0.02134 | | | 0.06
Pair | 4.1961 | 4.3786 | 4.5496 | 7.6 | 18.68
Bond | 0.00057323 | 0.00065148 | 0.00075008 | 0.0 | 0.00
Kspace | 11.998 | 12.168 | 12.35 | 4.6 | 51.91
Neigh | 0.064936 | 0.064967 | 0.064998 | 0.0 | 0.28
Comm | 0.089828 | 0.091261 | 0.092661 | 0.4 | 0.39
Output | 0.00082362 | 0.0010221 | 0.0015517 | 1.0 | 0.00
Modify | 6.7177 | 6.7195 | 6.7212 | 0.1 | 28.67
Other | | 0.01535 | | | 0.07
Nlocal: 944 ave 948 max 940 min
Histogram: 1 0 0 1 0 0 1 0 0 1
@ -168,4 +154,4 @@ Ave neighs/atom = 456.98835
Ave special neighs/atom = 0.50847458
Neighbor list builds = 6
Dangerous builds = 0
Total wall time: 0:01:01
Total wall time: 0:00:44

View File

@ -22,8 +22,8 @@
#include "compute.h"
#include "domain.h"
#include "electrode_accel_interface.h"
#include "electrode_matrix.h"
#include "electrode_math.h"
#include "electrode_matrix.h"
#include "electrode_vector.h"
#include "error.h"
#include "force.h"
@ -73,7 +73,7 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
ffield = false;
thermo_time = 0.;
scalar_flag = 1;
vector_flag = 1;
array_flag = 1;
top_group = 0;
intelflag = false;
tfflag = false;
@ -222,7 +222,8 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
"Cannot write elastance matrix if reading capacitance matrix "
"from file");
num_of_groups = static_cast<int>(groups.size());
size_vector = num_of_groups;
size_array_rows = num_of_groups;
size_array_cols = 2 + 2 * num_of_groups;
// check groups are consistent
int *mask = atom->mask;
@ -260,12 +261,6 @@ FixElectrodeConp::FixElectrodeConp(LAMMPS *lmp, int narg, char **arg) :
/* ---------------------------------------------------------------------- */
// 0 1 2 3 return 4 if successful
// fix_modify fxupdate set v [group] [var]
// fix_modify fxupdate set qsb [group] [var]
// 0 1 2 3 4 return 5 if successful
// fix_modify fxupdate set mc [group1] [group2] [var]
// fix_modify fxupdate set me [group1] [group2] [var]
int FixElectrodeConp::modify_param(int narg, char **arg)
{
if (strcmp(arg[0], "tf") == 0) {
@ -306,49 +301,6 @@ int FixElectrodeConp::modify_param(int narg, char **arg)
timer_flag = utils::logical(FLERR, arg[1], false, lmp);
return 2;
} else if (strcmp(arg[0], "set") == 0) {
if (strcmp(arg[1], "v") == 0 || strcmp(arg[1], "qsb") == 0) {
if (narg != 4)
error->all(FLERR,
fmt::format("Incorrect number of arguments for fix_modify set {}", arg[1]));
int outputgrp = groupnum_from_name(arg[2]);
int outputvar = input->variable->find(arg[3]);
if (outputvar < 0)
error->all(FLERR,
fmt::format("Variable {} for fix_modify electrode does not exist", arg[3]));
if (!input->variable->internalstyle(outputvar))
error->all(FLERR,
fmt::format("Variable {} for fix_modify electrode is invalid style", arg[3]));
if (strcmp(arg[1], "v") == 0)
setvars_types.push_back(V);
else
setvars_types.push_back(QSB);
setvars_groups.push_back(outputgrp);
setvars_vars.push_back(outputvar);
return 4;
} else if (strcmp(arg[1], "mc") == 0 || strcmp(arg[1], "me") == 0) {
if (narg != 5)
error->all(FLERR,
fmt::format("Incorrect number of arguments for fix_modify set {}", arg[1]));
int outputgrpi = groupnum_from_name(arg[2]);
int outputgrpj = groupnum_from_name(arg[3]);
int outputgrpij = outputgrpi * num_of_groups + outputgrpj;
int outputvar = input->variable->find(arg[4]);
if (outputvar < 0)
error->all(FLERR,
fmt::format("Variable {} for fix_modify electrode does not exist", arg[4]));
if (!input->variable->internalstyle(outputvar))
error->all(FLERR,
fmt::format("Variable {} for fix_modify electrode is invalid style", arg[4]));
if (strcmp(arg[1], "mc") == 0)
setvars_types.push_back(MC);
else
setvars_types.push_back(ME);
setvars_groups.push_back(outputgrpij);
setvars_vars.push_back(outputvar);
return 5;
} else
error->all(FLERR, "Invalid set option for fix_modify electrode");
} else
error->all(FLERR, "Invalid argument for fix_modify electrode");
return 0;
@ -452,10 +404,6 @@ void FixElectrodeConp::setup_post_neighbor()
int *mask = atom->mask;
tagint *tag = atom->tag;
// setvars asserts:
assert(setvars_groups.size() == setvars_vars.size());
assert(setvars_groups.size() == setvars_types.size());
// if Thomas-Fermi, make sure all electrode atoms have parameters
if (tfflag) {
int unset_tf = 0;
@ -539,8 +487,8 @@ void FixElectrodeConp::setup_post_neighbor()
compute_macro_matrices();
MPI_Barrier(world);
if (timer_flag && (comm->me == 0))
utils::logmesg(lmp,
fmt::format("SD-vector and macro matrices time: {:.4g} s\n", MPI_Wtime() - start));
utils::logmesg(
lmp, fmt::format("SD-vector and macro matrices time: {:.4g} s\n", MPI_Wtime() - start));
// initial charges and b vector
update_charges();
@ -789,9 +737,7 @@ void FixElectrodeConp::update_charges()
gather_elevec(charge_iele);
MPI_Allreduce(MPI_IN_PLACE, &sb_charges.front(), num_of_groups, MPI_DOUBLE, MPI_SUM, world);
update_setvars(QSB);
update_psi(); // use for equal-style and conq
update_setvars(V); // push psi into 'fix_modify set' vars
update_psi(); // use for equal-style and conq
for (int g = 0; g < num_of_groups; g++)
for (int j = 0; j < ngroup; j++) { charge_iele[j] += sd_vectors[g][j] * group_psi[g]; }
@ -820,34 +766,6 @@ void FixElectrodeConp::update_psi()
/* ---------------------------------------------------------------------- */
void FixElectrodeConp::update_setvars(int vtype)
{
int num_of_vars = static_cast<int>(setvars_groups.size());
for (int v = 0; v < num_of_vars; v++) {
if (vtype != setvars_types[v]) continue;
switch (vtype) {
case V:
input->variable->internal_set(setvars_vars[v], group_psi[setvars_groups[v]]);
break;
case QSB:
input->variable->internal_set(setvars_vars[v], sb_charges[setvars_groups[v]]);
break;
case MC: {
int groupi = setvars_groups[v] / num_of_groups;
int groupj = setvars_groups[v] % num_of_groups;
input->variable->internal_set(setvars_vars[v], macro_capacitance[groupi][groupj]);
} break;
case ME: {
int groupi = setvars_groups[v] / num_of_groups;
int groupj = setvars_groups[v] % num_of_groups;
input->variable->internal_set(setvars_vars[v], macro_elastance[groupi][groupj]);
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixElectrodeConp::compute_macro_matrices()
{
macro_capacitance =
@ -902,9 +820,6 @@ void FixElectrodeConp::compute_macro_matrices()
}
}
}
update_setvars(MC);
update_setvars(ME);
}
/* ---------------------------------------------------------------------- */
@ -915,9 +830,22 @@ double FixElectrodeConp::compute_scalar()
}
/* ---------------------------------------------------------------------- */
double FixElectrodeConp::compute_vector(int i)
double FixElectrodeConp::compute_array(int i, int j)
{
return group_psi[i];
if (i < 0 || i >= num_of_groups) error->all(FLERR, "invalid fix electrode array row reference");
if (j < 0)
error->all(FLERR, "invalid fix electrode array column reference");
else if (j == 0)
return group_psi[i];
else if (j == 1)
return sb_charges[i];
else if (j <= num_of_groups + 1)
return macro_capacitance[i][j - 2];
else if (j <= 2 * num_of_groups + 1)
return macro_elastance[i][j - num_of_groups - 2];
else
error->all(FLERR, "invalid fix electrode array column reference");
return 0.; // avoid -Wreturn-type warning
}
/* ---------------------------------------------------------------------- */

View File

@ -46,7 +46,7 @@ class FixElectrodeConp : public Fix {
void pre_force(int) override;
void pre_reverse(int, int) override;
double compute_scalar() override;
double compute_vector(int) override;
double compute_array(int, int) override;
int modify_param(int, char **) override;
int modify_param(const std::string &);
void init() override;
@ -104,8 +104,6 @@ class FixElectrodeConp : public Fix {
void read_from_file(std::string input_file, double **, const std::string &);
void compute_sd_vectors();
void compute_sd_vectors_ffield();
std::vector<int> setvars_types, setvars_groups, setvars_vars;
void update_setvars(int);
int groupnum_from_name(char *);
double evscale;
class Pair *pair;