update fix charge/regulation input example to follow LAMMPS conventions closer
This commit is contained in:
Axel Kohlmeyer
@ -1,8 +1,36 @@
|
||||
This directory has two input scripts that illustrates how to use fix charge_regulation in LAMMPS to perform coarse-grained molecular dynamics (MD) simulations with incorporation of charge regulation effects. The charge regulation is implemented via Monte Carlo (MC) sampling following the reaction ensemble MC approach, producing a MC/MD hybrid tool for modeling charge regulation in solvated systems.
|
||||
This directory has two input scripts that illustrates how to use fix
|
||||
charge_regulation in LAMMPS to perform coarse-grained molecular dynamics
|
||||
(MD) simulations with incorporation of charge regulation effects. The
|
||||
charge regulation is implemented via Monte Carlo (MC) sampling following
|
||||
the reaction ensemble MC approach, producing a MC/MD hybrid tool for
|
||||
modeling charge regulation in solvated systems.
|
||||
|
||||
The script `in_acid.chreg` sets up a simple weak acid electrolyte (pH=7,pKa=6,pI=3). Four different types of MC moves are implemented: acid protonation & de-protonation, and monovalent ion pair insertion and deletion. Note here we have grouped all free monovalent ions into a single type, a physically natural choice on the level of coarse-grained primitive electrolyte models, which increases the calculation performance but has no effects on thermodynamic observables. The variables such as pH, pKa, pI, and lb at the top of the input script can be adjusted to play with various simulation parameters. The cumulative MC attempted moves and cumulative number of accepted moves, as well as, current number of neutral and charged acid particles, neutral and charged base particles (in this example always 0), and the current number of free cations and anions in the system are printed in the `log_acid.lammps`.
|
||||
The script `in.chreg-acid` sets up a simple weak acid electrolyte
|
||||
(pH=7,pKa=6,pI=3). Four different types of MC moves are implemented:
|
||||
acid protonation & de-protonation, and monovalent ion pair insertion and
|
||||
deletion. Note here we have grouped all free monovalent ions into a
|
||||
single type, a physically natural choice on the level of coarse-grained
|
||||
primitive electrolyte models, which increases the calculation
|
||||
performance but has no effects on thermodynamic observables. The
|
||||
variables such as pH, pKa, pI, and lb at the top of the input script can
|
||||
be adjusted to play with various simulation parameters. The cumulative
|
||||
MC attempted moves and cumulative number of accepted moves, as well as,
|
||||
current number of neutral and charged acid particles, neutral and
|
||||
charged base particles (in this example always 0), and the current
|
||||
number of free cations and anions in the system are printed in the
|
||||
output.
|
||||
|
||||
The script `in_polymer.chreg` sets up a weak polyelectrolyte chain of N=80 beads. Each bead is a weak acid with pKa=5 and solution has pH=7 and monovalent salt chemical potential pI=3. In this example, we choose to treat salt ions, protons, and hydroxyl ions separately, which results in 5 types of MC moves: acid [type 1] protonation & de-protonation (with protons [type 4] insertion & deletion), acid [type 1] protonation & de-protonation (with salt cation [type 2] insertion & deletion), water self-ionization (insertion and deletion of proton [type4] and hydroxyl ion [type 5] pair), insertion and deletion of monovalent salt pair [type 2 and type 3] , insertion and deletion
|
||||
Of a proton [type4] and salt anion [type 3].
|
||||
The current number of neutral and charged acid particles, the current number of free salt cations and anions, and the current number of protons and hydroxyl ions are printed in the `log_polymer.lammps`.
|
||||
The script `in.chreg-polymer` sets up a weak poly-electrolyte chain of
|
||||
N=80 beads. Each bead is a weak acid with pKa=5 and solution has pH=7
|
||||
and monovalent salt chemical potential pI=3. In this example, we choose
|
||||
to treat salt ions, protons, and hydroxyl ions separately, which results
|
||||
in 5 types of MC moves: acid [type 1] protonation & de-protonation (with
|
||||
protons [type 4] insertion & deletion), acid [type 1] protonation &
|
||||
de-protonation (with salt cation [type 2] insertion & deletion), water
|
||||
self-ionization (insertion and deletion of proton [type4] and hydroxyl
|
||||
ion [type 5] pair), insertion and deletion of monovalent salt pair [type
|
||||
2 and type 3] , insertion and deletion of a proton [type4] and salt
|
||||
anion [type 3]. The current number of neutral and charged acid
|
||||
particles, the current number of free salt cations and anions, and the
|
||||
current number of protons and hydroxyl ions are printed in the output.
|
||||
|
||||
|
||||
@ -1,12 +1,9 @@
|
||||
# Charge regulation lammps for simple weak electrolyte
|
||||
|
||||
units lj
|
||||
atom_style charge
|
||||
dimension 3
|
||||
boundary p p p
|
||||
processors * * *
|
||||
neighbor 3.0 bin
|
||||
read_data data_acid.chreg
|
||||
units lj
|
||||
atom_style charge
|
||||
neighbor 3.0 bin
|
||||
read_data data.chreg-acid
|
||||
|
||||
variable cut_long equal 12.5
|
||||
variable nevery equal 100
|
||||
@ -16,26 +13,24 @@ variable pKa equal 6.0
|
||||
variable pIm equal 3.0
|
||||
variable pIp equal 3.0
|
||||
|
||||
variable cut_lj equal 2^(1.0/6.0)
|
||||
variable lunit_nm equal 0.72 # in the units of nm
|
||||
velocity all create 1.0 8008 loop geom
|
||||
variable cut_lj equal 2^(1.0/6.0)
|
||||
variable lunit_nm equal 0.72 # in the units of nm
|
||||
velocity all create 1.0 8008 loop geom
|
||||
|
||||
pair_style lj/cut/coul/long ${cut_lj} ${cut_long}
|
||||
pair_coeff * * 1.0 1.0 ${cut_lj} # charges
|
||||
kspace_style ewald 1.0e-3
|
||||
dielectric 1.0
|
||||
pair_modify shift yes
|
||||
pair_coeff * * 1.0 1.0
|
||||
kspace_style ewald 1.0e-3
|
||||
pair_modify shift yes
|
||||
|
||||
######### VERLET INTEGRATION WITH LANGEVIN THERMOSTAT ###########
|
||||
fix fnve all nve
|
||||
fix fnve all nve
|
||||
compute dtemp all temp
|
||||
compute_modify dtemp dynamic yes
|
||||
fix fT all langevin 1.0 1.0 1.0 123
|
||||
fix_modify fT temp dtemp
|
||||
|
||||
fix chareg all charge/regulation 2 3 acid_type 1 pH ${pH} pKa ${pKa} pIp ${pIp} pIm ${pIm} lunit_nm ${lunit_nm} nevery ${nevery} nmc ${nmc} seed 2345 tempfixid fT
|
||||
thermo 100
|
||||
thermo 100
|
||||
thermo_style custom step pe c_dtemp f_chareg[1] f_chareg[2] f_chareg[3] f_chareg[4] f_chareg[5] f_chareg[6] f_chareg[7] f_chareg[8]
|
||||
log log_acid.lammps
|
||||
timestep 0.005
|
||||
run 20000
|
||||
timestep 0.005
|
||||
run 2000
|
||||
33
examples/USER/misc/charge_regulation/in.chreg-polymer
Normal file
33
examples/USER/misc/charge_regulation/in.chreg-polymer
Normal file
@ -0,0 +1,33 @@
|
||||
# Charge regulation lammps for a polymer chain
|
||||
units lj
|
||||
atom_style full
|
||||
neighbor 3.0 bin
|
||||
read_data data.chreg-polymer
|
||||
|
||||
bond_style harmonic
|
||||
bond_coeff 1 100 1.122462 # K R0
|
||||
velocity all create 1.0 8008 loop geom
|
||||
|
||||
pair_style lj/cut/coul/long 1.122462 20
|
||||
pair_coeff * * 1.0 1.0 1.122462 # charges
|
||||
kspace_style pppm 1.0e-3
|
||||
pair_modify shift yes
|
||||
dielectric 1.0
|
||||
|
||||
######### VERLET INTEGRATION WITH LANGEVIN THERMOSTAT ###########
|
||||
fix fnve all nve
|
||||
compute dtemp all temp
|
||||
compute_modify dtemp dynamic yes
|
||||
fix fT all langevin 1.0 1.0 1.0 123
|
||||
fix_modify fT temp dtemp
|
||||
|
||||
fix chareg1 all charge/regulation 2 3 acid_type 1 pH 7.0 pKa 6.5 pIp 3.0 pIm 3.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg2 all charge/regulation 4 5 acid_type 1 pH 7.0 pKa 6.5 pIp 7.0 pIm 7.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg3 all charge/regulation 4 3 pIp 7.0 pIm 3.0 temp 1.0 nmc 20 seed 2345
|
||||
|
||||
thermo 100
|
||||
# print: step, potential energy, temperature, neutral acids, charged acids, salt cations, salt anions, H+ ions, OH- ions
|
||||
thermo_style custom step pe c_dtemp f_chareg1[3] f_chareg1[4] f_chareg1[7] f_chareg1[8] f_chareg2[7] f_chareg2[8]
|
||||
|
||||
timestep 0.005
|
||||
run 2000
|
||||
@ -1,37 +0,0 @@
|
||||
# Charge regulation lammps for a polymer chain
|
||||
units lj
|
||||
atom_style full
|
||||
dimension 3
|
||||
boundary p p p
|
||||
processors * * *
|
||||
neighbor 3.0 bin
|
||||
read_data data_polymer.chreg
|
||||
|
||||
bond_style harmonic
|
||||
bond_coeff 1 100 1.122462 # K R0
|
||||
velocity all create 1.0 8008 loop geom
|
||||
|
||||
pair_style lj/cut/coul/long 1.122462 20
|
||||
pair_coeff * * 1.0 1.0 1.122462 # charges
|
||||
kspace_style pppm 1.0e-3
|
||||
pair_modify shift yes
|
||||
dielectric 1.0
|
||||
|
||||
######### VERLET INTEGRATION WITH LANGEVIN THERMOSTAT ###########
|
||||
fix fnve all nve
|
||||
compute dtemp all temp
|
||||
compute_modify dtemp dynamic yes
|
||||
fix fT all langevin 1.0 1.0 1.0 123
|
||||
fix_modify fT temp dtemp
|
||||
|
||||
fix chareg1 all charge/regulation 2 3 acid_type 1 pH 7.0 pKa 6.5 pIp 3.0 pIm 3.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg2 all charge/regulation 4 5 acid_type 1 pH 7.0 pKa 6.5 pIp 7.0 pIm 7.0 temp 1.0 nmc 40 seed 2345
|
||||
fix chareg3 all charge/regulation 4 3 pIp 7.0 pIm 3.0 temp 1.0 nmc 20 seed 2345
|
||||
|
||||
thermo 100
|
||||
# print: step, potential energy, temperature, neutral acids, charged acids, salt cations, salt anions, H+ ions, OH- ions
|
||||
thermo_style custom step pe c_dtemp f_chareg1[3] f_chareg1[4] f_chareg1[7] f_chareg1[8] f_chareg2[7] f_chareg2[8]
|
||||
log log_polymer.lammps
|
||||
|
||||
timestep 0.005
|
||||
run 20000
|
||||
@ -1,163 +0,0 @@
|
||||
LAMMPS (24 Dec 2020)
|
||||
Reading data file ...
|
||||
orthogonal box = (-25.000000 -25.000000 -25.000000) to (25.000000 25.000000 25.000000)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
219 atoms
|
||||
read_data CPU = 0.001 seconds
|
||||
Ewald initialization ...
|
||||
using 12-bit tables for long-range coulomb (../kspace.cpp:339)
|
||||
G vector (1/distance) = 0.14221027
|
||||
estimated absolute RMS force accuracy = 0.0010128126
|
||||
estimated relative force accuracy = 0.0010128126
|
||||
KSpace vectors: actual max1d max3d = 257 5 665
|
||||
kxmax kymax kzmax = 5 5 5
|
||||
0 atoms in group Fix_CR:exclusion_group:chareg
|
||||
WARNING: Neighbor exclusions used with KSpace solver may give inconsistent Coulombic energies (../neighbor.cpp:486)
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 10 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 15.5
|
||||
ghost atom cutoff = 15.5
|
||||
binsize = 7.75, bins = 7 7 7
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cut/coul/long, perpetual
|
||||
attributes: half, newton on
|
||||
pair build: half/bin/atomonly/newton
|
||||
stencil: half/bin/3d/newton
|
||||
bin: standard
|
||||
Setting up Verlet run ...
|
||||
Unit style : lj
|
||||
Current step : 0
|
||||
Time step : 0.005
|
||||
Per MPI rank memory allocation (min/avg/max) = 11.91 | 11.91 | 11.91 Mbytes
|
||||
Step PotEng c_dtemp f_chareg[1] f_chareg[2] f_chareg[3] f_chareg[4] f_chareg[5] f_chareg[6] f_chareg[7] f_chareg[8]
|
||||
0 -0.054228269 1 0 0 1 99 0 0 109 10
|
||||
100 -0.058672881 0.99159291 100 71 16 84 0 0 92 8
|
||||
200 -0.05399313 0.92006523 200 154 26 74 0 0 85 11
|
||||
300 -0.047035343 0.92728143 300 240 22 78 0 0 85 7
|
||||
400 -0.049597809 1.0617937 400 319 16 84 0 0 92 8
|
||||
500 -0.052409835 1.0006148 500 404 12 88 0 0 97 9
|
||||
600 -0.056012172 0.98059344 600 481 15 85 0 0 92 7
|
||||
700 -0.053639989 0.99572709 700 561 16 84 0 0 94 10
|
||||
800 -0.060026132 0.95764632 800 639 22 78 0 0 84 6
|
||||
900 -0.050785422 0.98399084 900 719 28 72 0 0 82 10
|
||||
1000 -0.062294743 0.97200068 1000 797 26 74 0 0 82 8
|
||||
1100 -0.051269402 1.0064376 1100 877 25 75 0 0 84 9
|
||||
1200 -0.077744839 1.0159098 1200 955 23 77 0 0 88 11
|
||||
1300 -0.084889696 1.1230485 1300 1037 20 80 0 0 90 10
|
||||
1400 -0.059361445 0.96735845 1400 1120 18 82 0 0 93 11
|
||||
1500 -0.052926174 0.95579188 1500 1199 24 76 0 0 86 10
|
||||
1600 -0.052376649 0.99376378 1600 1284 22 78 0 0 89 11
|
||||
1700 -0.052480188 0.96085964 1700 1361 27 73 0 0 84 11
|
||||
1800 -0.065884306 0.96747971 1800 1441 21 79 0 0 84 5
|
||||
1900 -0.054315859 0.95873145 1900 1521 23 77 0 0 84 7
|
||||
2000 -0.037161802 0.93562039 2000 1604 27 73 0 0 79 6
|
||||
2100 -0.034977265 0.97177103 2100 1684 26 74 0 0 85 11
|
||||
2200 -0.047434868 0.97897613 2200 1762 17 83 0 0 93 10
|
||||
2300 -0.047392634 0.96570672 2300 1837 18 82 0 0 92 10
|
||||
2400 -0.044879306 0.98620033 2400 1910 19 81 0 0 89 8
|
||||
2500 -0.069690496 1.0690505 2500 1988 16 84 0 0 91 7
|
||||
2600 -0.081588407 0.97711054 2600 2067 17 83 0 0 92 9
|
||||
2700 -0.06341681 1.0386711 2700 2146 20 80 0 0 85 5
|
||||
2800 -0.045290012 1.0402055 2800 2230 18 82 0 0 91 9
|
||||
2900 -0.046875012 1.0609775 2900 2317 24 76 0 0 86 10
|
||||
3000 -0.031258722 0.93861202 3000 2400 24 76 0 0 85 9
|
||||
3100 -0.04673342 0.90800583 3100 2485 25 75 0 0 90 15
|
||||
3200 -0.054354227 0.94493881 3200 2567 16 84 0 0 94 10
|
||||
3300 -0.053647746 0.92321446 3300 2641 17 83 0 0 94 11
|
||||
3400 -0.031751732 0.93735127 3400 2725 22 78 0 0 92 14
|
||||
3500 -0.053806113 0.98798136 3500 2795 28 72 0 0 84 12
|
||||
3600 -0.040751349 0.84291639 3600 2873 28 72 0 0 84 12
|
||||
3700 -0.051747138 1.072448 3700 2951 24 76 0 0 92 16
|
||||
3800 -0.043420594 1.0076309 3800 3030 26 74 0 0 79 5
|
||||
3900 -0.050845848 1.0250023 3900 3103 29 71 0 0 76 5
|
||||
4000 -0.039837847 1.064111 4000 3182 29 71 0 0 83 12
|
||||
4100 -0.045638995 1.1249685 4100 3262 28 72 0 0 81 9
|
||||
4200 -0.047956491 0.92255907 4200 3348 26 74 0 0 87 13
|
||||
4300 -0.054052822 1.006239 4300 3423 19 81 0 0 90 9
|
||||
4400 -0.053148034 1.0028887 4400 3506 24 76 0 0 83 7
|
||||
4500 -0.062132076 1.0317847 4500 3587 23 77 0 0 82 5
|
||||
4600 -0.04616043 0.99066453 4600 3673 28 72 0 0 82 10
|
||||
4700 -0.066990889 1.0242064 4700 3755 18 82 0 0 90 8
|
||||
4800 -0.0564736 0.91765628 4800 3832 22 78 0 0 91 13
|
||||
4900 -0.052301294 0.95348659 4900 3912 28 72 0 0 81 9
|
||||
5000 -0.062630677 1.0336579 5000 3989 18 82 0 0 92 10
|
||||
5100 -0.053645908 1.0314613 5100 4062 20 80 0 0 85 5
|
||||
5200 -0.062189568 1.0504732 5200 4133 23 77 0 0 84 7
|
||||
5300 -0.049092746 1.0310832 5300 4217 24 76 0 0 82 6
|
||||
5400 -0.051859257 0.99600428 5400 4299 27 73 0 0 80 7
|
||||
5500 -0.065540815 0.98012128 5500 4381 19 81 0 0 92 11
|
||||
5600 -0.071018582 0.9252814 5600 4455 23 77 0 0 84 7
|
||||
5700 -0.066954185 1.0325214 5700 4535 26 74 0 0 82 8
|
||||
5800 -0.064847249 1.0313536 5800 4617 21 79 0 0 90 11
|
||||
5900 -0.063173056 0.95455853 5900 4703 18 82 0 0 92 10
|
||||
6000 -0.064807837 0.97182222 6000 4790 21 79 0 0 91 12
|
||||
6100 -0.07067683 0.91775921 6100 4875 16 84 0 0 94 10
|
||||
6200 -0.071400842 1.0162225 6200 4952 17 83 0 0 87 4
|
||||
6300 -0.078479449 1.0106706 6300 5033 21 79 0 0 84 5
|
||||
6400 -0.083167651 1.0246584 6400 5111 18 82 0 0 90 8
|
||||
6500 -0.092611537 1.0766467 6500 5195 20 80 0 0 88 8
|
||||
6600 -0.096710071 1.0246648 6600 5274 15 85 0 0 91 6
|
||||
6700 -0.073399857 0.94939392 6700 5351 17 83 0 0 92 9
|
||||
6800 -0.062479375 0.9393967 6800 5434 18 82 0 0 93 11
|
||||
6900 -0.065391043 0.93475954 6900 5514 22 78 0 0 89 11
|
||||
7000 -0.045655499 0.98688239 7000 5601 21 79 0 0 90 11
|
||||
7100 -0.061186309 1.0309063 7100 5684 22 78 0 0 87 9
|
||||
7200 -0.074737514 1.0516593 7200 5769 25 75 0 0 80 5
|
||||
7300 -0.075228979 1.0167704 7300 5847 21 79 0 0 86 7
|
||||
7400 -0.06660147 1.0947107 7400 5930 25 75 0 0 87 12
|
||||
7500 -0.071915247 1.10542 7500 6009 24 76 0 0 84 8
|
||||
7600 -0.029974303 0.99202697 7600 6095 28 72 0 0 80 8
|
||||
7700 -0.029024004 1.0390995 7700 6182 28 72 0 0 83 11
|
||||
7800 -0.056250746 0.96393961 7800 6260 18 82 0 0 91 9
|
||||
7900 -0.072178944 0.9833919 7900 6339 21 79 0 0 86 7
|
||||
8000 -0.070377248 1.021342 8000 6419 23 77 0 0 88 11
|
||||
8100 -0.07753283 0.96194312 8100 6493 26 74 0 0 86 12
|
||||
8200 -0.075617309 0.9715344 8200 6576 25 75 0 0 89 14
|
||||
8300 -0.077480013 0.99106705 8300 6662 20 80 0 0 91 11
|
||||
8400 -0.079901928 0.89855112 8400 6744 23 77 0 0 91 14
|
||||
8500 -0.069192745 1.0606791 8500 6822 19 81 0 0 91 10
|
||||
8600 -0.058657202 1.1270217 8600 6898 15 85 0 0 95 10
|
||||
8700 -0.044985397 1.0367419 8700 6979 27 73 0 0 84 11
|
||||
8800 -0.064266376 0.91557003 8800 7060 21 79 0 0 89 10
|
||||
8900 -0.068680533 0.95963643 8900 7133 19 81 0 0 88 7
|
||||
9000 -0.051736144 1.0398547 9000 7213 13 87 0 0 94 7
|
||||
9100 -0.058381436 0.98047663 9100 7290 17 83 0 0 89 6
|
||||
9200 -0.05531014 0.92541631 9200 7368 22 78 0 0 85 7
|
||||
9300 -0.04386907 0.9339658 9300 7454 22 78 0 0 89 11
|
||||
9400 -0.052293168 0.94314034 9400 7539 22 78 0 0 86 8
|
||||
9500 -0.050362046 1.0489028 9500 7617 18 82 0 0 90 8
|
||||
9600 -0.054272227 1.0879161 9600 7697 11 89 0 0 96 7
|
||||
9700 -0.042514179 1.0051505 9700 7771 22 78 0 0 84 6
|
||||
9800 -0.045365018 0.95363669 9800 7850 25 75 0 0 77 2
|
||||
9900 -0.064287274 0.91994667 9900 7928 27 73 0 0 81 8
|
||||
10000 -0.05689162 0.99963208 10000 8011 22 78 0 0 84 6
|
||||
Loop time of 19.4452 on 1 procs for 10000 steps with 190 atoms
|
||||
|
||||
Performance: 222162.510 tau/day, 514.265 timesteps/s
|
||||
99.9% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.68864 | 0.68864 | 0.68864 | 0.0 | 3.54
|
||||
Kspace | 6.7893 | 6.7893 | 6.7893 | 0.0 | 34.92
|
||||
Neigh | 0.060782 | 0.060782 | 0.060782 | 0.0 | 0.31
|
||||
Comm | 0.047035 | 0.047035 | 0.047035 | 0.0 | 0.24
|
||||
Output | 0.0027227 | 0.0027227 | 0.0027227 | 0.0 | 0.01
|
||||
Modify | 11.838 | 11.838 | 11.838 | 0.0 | 60.88
|
||||
Other | | 0.01878 | | | 0.10
|
||||
|
||||
Nlocal: 190.000 ave 190 max 190 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 592.000 ave 592 max 592 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 2194.00 ave 2194 max 2194 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2194
|
||||
Ave neighs/atom = 11.547368
|
||||
Neighbor list builds = 10287
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:19
|
||||
Reference in New Issue
Block a user