Dear LAMMPS users
For a testing & reference purpose, I have to run a test simulation of one single Na+ in a vacuum box of various box sizes to compute the free energy of charging/discharging the ion under the PBC conditions.
I know that generally Ewald summation and related methods can deal with charge non-neutrality as if in a neutralizing background and hence this should be case in LAMMPS as well (but there is warning)
I am using the latest 28th June version and I have the following input file and simple initial configuration (actually I did not make it run but just output the energy with the ion centered in the cubic box in a few steps).
For this setting, the total pe (which is also the elong output as the ionic pair across the box is beyond cutoff 12 Armstrong from each other) is ~ -18 kcal/mol-1, which means that the free energy of charging the ion from neutral to charge +1 under such PBC (L = 25.7Armstrong) should be around ~ - 18kcal mol-1 (indeed as I have computed).
This number seems too large for a Ewald PBC finite-size correction to the charging free energy (compare with the free energy of charging an ion in a real vacuum should be = 0) as various reference states that it should be around ~ 1kcal mol-1, or at least not a large number
I might be mistaken for a lot of concepts and information, but I have the following initial configuration and input file.
Many thanks!
- config.dat
LAMMPS Description
1 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
1 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
-12.7 12.7 xlo xhi
-12.7 12.7 ylo yhi
-12.7 12.7 zlo zhi
Masses
1 22.9898
Atoms
1 1 1 1.00000 0.00 0.00 0.00 0 0 0
2）input.dat
#restart from the only initiation binary
dimension 3
units real
atom_style full
boundary p p p
read_data config-naalone.dat
#replicate 2 2 2
group all type 1
group Na type 1
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 4000 #exclude molecule all
set group Na charge 1.000
pair_style lj/cut/coul/long 12.0 12.0
pair_coeff 1 1 0.0 0.0
kspace_style pppm 1.0e-4
#kspace_modify gewald 0.26
pair_modify shift no
pair_modify tail yes
#velocity all create 298 4271 dist gaussian mom yes rot yes
#fix 2 all nvt temp 298.0 298.0 100
#fix 3 all recenter INIT INIT INIT
restart 1000 old_config_1.dat old_config_2.dat
compute 14 all com
thermo 1000
thermo_style custom step temp pe density evdwl etail ecoul elong press vol c_14[1] c_14[2] c_14[3] lx ly lz emol ke
dump id all xyz 1000 dump_init.xyz
dump_modify id element Na
timestep 1.0
run 5
write_data data_restart_init
3）output
LAMMPS (28 Jun 2014)
#restart from the only initiation binary
dimension 3
units real
atom_style full
boundary p p p
read_data config-naalone.dat
orthogonal box = (-12.7 -12.7 -12.7) to (12.7 12.7 12.7)
1 by 1 by 1 MPI processor grid
reading atoms …
1 atoms
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
#replicate 2 2 2
group all type 1
1 atoms in group all
group Na type 1
1 atoms in group Na
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes one 4000 #exclude molecule all
set group Na charge 1.000
1 settings made for charge
pair_style lj/cut/coul/long 12.0 12.0
pair_coeff 1 1 0.0 0.0
kspace_style pppm 1.0e-4
#kspace_modify gewald 0.26
pair_modify shift no
pair_modify tail yes
#velocity all create 298 4271 dist gaussian mom yes rot yes
#fix 2 all nvt temp 298.0 298.0 100
#fix 3 all recenter INIT INIT INIT
restart 1000 old_config_1.dat old_config_2.dat
compute 14 all com
thermo 1000
thermo_style custom step temp pe density evdwl etail ecoul elong press vol c_14[1] c_14[2] c_14[3] lx ly lz emol ke
dump id all xyz 1000 dump_init.xyz
dump_modify id element Na
timestep 1.0
run 5
WARNING: No fixes defined, atoms won’t move (…/verlet.cpp:55)
PPPM initialization …
WARNING: System is not charge neutral, net charge = 1 (…/kspace.cpp:278)
G vector (1/distance) = 0.196752
grid = 8 8 8
stencil order = 5
estimated absolute RMS force accuracy = 0.00803509
estimated relative force accuracy = 2.41974e-05
using double precision FFTs
3d grid and FFT values/proc = 2197 512
Memory usage per processor = 6.31464 Mbytes
Step Temp PotEng Density E_vdwl E_tail E_coul E_long Press Volume 14[1] 14[2] 14[3] Lx Ly Lz E_mol KinEng
0 0 -18.486529 0.0023296092 0 0 0 -18.486529 -22.863308 16387.064 0 0 0 25.4
25.4 25.4 0 0
5 0 -18.486529 0.0023296092 0 0 0 -18.486529 -22.863308 16387.064 0 0 0 25.4
25.4 25.4 0 0
Loop time of 0.000310898 on 1 procs for 5 steps with 1 atoms
Pair time () = 1.19209e-06 (0.383436)
Bond time () = 0 (0)
Kspce time () = 0.000285625 (91.8712)
Neigh time () = 0 (0)
Comm time () = 7.15256e-06 (2.30061)
Outpt time () = 1.50204e-05 (4.83129)
Other time (%) = 1.90735e-06 (0.613497)
Nlocal: 1 ave 1 max 1 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 26 ave 26 max 26 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 0
Ave neighs/atom = 0
Ave special neighs/atom = 0
Neighbor list builds = 0
Dangerous builds = 0
write_data data_restart_init
PPPM initialization …
WARNING: System is not charge neutral, net charge = 1 (…/kspace.cpp:278)
G vector (1/distance) = 0.196752
grid = 8 8 8
stencil order = 5
estimated absolute RMS force accuracy = 0.00803509
estimated relative force accuracy = 2.41974e-05
using double precision FFTs
3d grid and FFT values/proc = 2197 512
Best Regards
Lunna