84 lines
2.6 KiB
Plaintext
84 lines
2.6 KiB
Plaintext
# Two ions, a cation and an anion, confined between two interfaces: epsilon1 | epsilon2 | epsilon1
|
|
# The interface normal vectors should be consistent with ed, pointing from region with epsilon1 to that with epsilon2
|
|
# bottom interface: n = (0, 0, 1)
|
|
# top interface: n = (0, 0, -1)
|
|
# so that ed's are the same for both interfaces
|
|
|
|
variable epsilon1 index 20
|
|
variable epsilon2 index 8
|
|
|
|
variable data index data.confined
|
|
|
|
newton off
|
|
units lj
|
|
atom_style dielectric
|
|
atom_modify map array
|
|
dimension 3
|
|
boundary p p f
|
|
|
|
variable method index gmres # gmres = BEM/GMRES
|
|
# icc = BEM/ICC*
|
|
# dof = Direct optimization of the functional
|
|
# none
|
|
|
|
variable ed equal "v_epsilon2 - v_epsilon1"
|
|
variable em equal "(v_epsilon2 + v_epsilon1)/2"
|
|
variable epsilon equal 1.0 # epsilon at the patch, not used for now
|
|
variable area equal 0.866 # patch area, same as in the data file
|
|
|
|
read_data ${data}
|
|
|
|
group interface type 1
|
|
group ions type 2 3
|
|
|
|
group cations type 2
|
|
group anions type 3
|
|
|
|
# 1.0 = q * epsilon2 = qreal for cations
|
|
# -1.0 = q * epsilon2 = qreal for anions
|
|
variable qscale equal "1.0 / v_epsilon2"
|
|
set group cations charge ${qscale}
|
|
variable qscale equal "-1.0 / v_epsilon2"
|
|
set group anions charge ${qscale}
|
|
|
|
pair_style lj/cut/coul/long/dielectric 1.122 10.0
|
|
pair_coeff * * 1.0 1.0
|
|
pair_coeff 1 1 0.0 1.0
|
|
|
|
kspace_style pppm/dielectric 0.0001
|
|
kspace_modify slab 3.0
|
|
|
|
neigh_modify every 1 delay 0 check yes one 5000
|
|
|
|
#compute ef all efield/atom
|
|
dump 1 all custom 100 all.dump id mol type q x y z #fx fy fz c_ef[1] c_ef[2] c_ef[3]
|
|
dump 2 interface custom 100 interface.dump id mol type q x y z #fx fy fz c_ef[1] c_ef[2] c_ef[3]
|
|
dump_modify 1 sort id
|
|
|
|
dump 3 ions custom 100 ions.dump id mol type q x y z fx fy fz #c_ef[1] c_ef[2] c_ef[3]
|
|
|
|
fix 1 ions nve
|
|
|
|
if "${method} == gmres" then &
|
|
"fix 3 interface polarize/bem/gmres 1 1.0e-4" &
|
|
"fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &
|
|
elif "${method} == icc"&
|
|
"fix 3 interface polarize/bem/icc 1 1.0e-4 itr_max 50" &
|
|
"fix_modify 3 itr_max 50 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &
|
|
elif "${method} == dof" &
|
|
"fix 3 interface polarize/functional 1 0.001" &
|
|
"fix_modify 3 dielectrics ${ed} ${em} ${epsilon} ${area} NULL" &
|
|
else &
|
|
"print 'Unsupported polarization solver' "
|
|
|
|
thermo 1000
|
|
thermo_style custom step evdwl ecoul elong epair #f_3
|
|
thermo_modify flush yes
|
|
|
|
run 0
|
|
|
|
|
|
|
|
|
|
|