I am attempting to verify the performance of the fix efield command. I set up a simple system with a few TIP3P water molecules. I calculated the forces on each atom with no external field and with an applied external field of .1 V/angstrom in the direction of the z-axis. Then, I took the difference in z component of the forces on each atom (with and without an external field) and divided by the charge of that atom. After carrying out the appropriate unit conversions this should yield a value of .1 V/angstrom, but I am getting a value of ~.26. Is this an issue with the implementation, or am I making a mistake? Lammps input files are shown below as well as a python script to calculate the applied field.
----------- structure.dat ---------------
LAMMPS data file via write_data, version 12 Dec 2018, timestep = 19240000
9 atoms
2 atom types
6 bonds
1 bond types
3 angles
1 angle types
0 impropers
0 improper types
1.1703503591503361e+00 4.7829649640855514e+01 xlo xhi
1.1703503591503361e+00 4.7829649640855514e+01 ylo yhi
1.1703503591503361e+00 4.7829649640855514e+01 zlo zhi
Masses
1 15.9994
2 1.00794
Pair Coeffs # lj/cut/coul/long
1 0.102 3.188
2 0 0
Bond Coeffs # harmonic
1 450 0.9572
Angle Coeffs # harmonic
1 55 104.52
Atoms # full
196 192 1 -8.3399999999999996e-01 4.6633003919069488e+01 1.4340843855491901e+01 3.1773569106055568e+00 -2 -7 3
194 192 2 4.1699999999999998e-01 4.7047863995991385e+01 1.4717487311993107e+01 3.9534127483204888e+00 -2 -7 3
195 192 2 4.1699999999999998e-01 4.7031741570866458e+01 1.4806469312634462e+01 2.4422153876993460e+00 -2 -7 3
193 191 1 -8.3399999999999996e-01 4.6145767669394068e+01 6.8469346709683601e+00 2.5403089509668018e+01 0 -3 -8
192 191 2 4.1699999999999998e-01 4.5833390713831839e+01 6.0885689626130448e+00 2.5896581080036292e+01 0 -3 -8
191 191 2 4.1699999999999998e-01 4.7091824283587300e+01 6.7171563113153105e+00 2.5337011524625833e+01 0 -3 -8
199 193 1 -8.3399999999999996e-01 3.9011683390379332e+01 2.5033409568055053e+01 4.5593393738722817e+01 -4 4 -1
197 193 2 4.1699999999999998e-01 3.9253185355361985e+01 2.5814683714860536e+01 4.6090906902813018e+01 -4 4 -1
198 193 2 4.1699999999999998e-01 3.9399947352470264e+01 2.5169708698958217e+01 4.4729156802080553e+01 -4 4 -1
Velocities
196 7.3121330182025605e-04 2.4659344353605678e-03 2.5949486192051716e-03
194 9.4401731689048012e-03 -2.4200211680022285e-04 -7.5648134923503513e-04
195 1.4529464690668983e-03 -3.6451415952267687e-03 -8.6719269555189096e-04
193 -1.2783819273667484e-03 -6.0529701689978309e-03 2.8400417261422060e-03
192 4.6824849880882976e-03 1.3295216559454132e-03 1.7941305077018601e-02
191 7.0095410126012837e-04 3.7938190301906277e-03 1.1758612330964082e-02
199 -5.2870289341674639e-03 8.9154982785344763e-04 2.2501048282937606e-03
197 5.9487474203909962e-03 -6.1718506498871584e-03 7.8904274107742811e-03
198 9.7344352511807928e-03 -5.5309951927488337e-03 7.9966492684395169e-03
Bonds
271 1 196 194
272 1 196 195
3061 1 193 191
3062 1 193 192
6381 1 199 197
6382 1 199 198
Angles
136 1 195 196 194
1531 1 192 193 191
3191 1 198 199 197
-------------- inputNoEfield ---------------------------
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
improper_style harmonic
pair_style lj/cut/coul/long 10.000000000000000
read_data structure.dat
run_style verlet
timestep 1.000000000000000
pair_coeff 1 1 0.102000000000000 3.188000000000000 10.000000000000000
pair_coeff 1 2 0.000000000000000 0.000000000000000 10.000000000000000
pair_coeff 2 2 0.000000000000000 0.000000000000000 10.000000000000000
dielectric 1.0
kspace_style pppm 0.000100000000000
mass 1 15.999400000000000
mass 2 1.007940000000000
bond_coeff 1 450.000000000000000 0.957200000000000
angle_coeff 1 55.000000000000000 104.519999999999996
fix mobileEnsStr all nvt temp 300.000000000000000 300.000000000000000 100.000000000000000
fix shakeFixB1A1GwaterGrpStr all shake 0.0001 100 0 b 1 a 1
timestep 1.0
dump dump1 all custom 1 noEfield.dump id mol q x y z fx fy fz
dump_modify dump1 sort id
run 0
------- inputWithEfield -----------------------
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
pair_style lj/cut/coul/long 10.000000000000000
read_data structure.dat
run_style verlet
timestep 1.000000000000000
pair_coeff 1 1 0.102000000000000 3.188000000000000 10.000000000000000
pair_coeff 1 2 0.000000000000000 0.000000000000000 10.000000000000000
pair_coeff 2 2 0.000000000000000 0.000000000000000 10.000000000000000
dielectric 1.0
kspace_style pppm 0.000100000000000
mass 1 15.999400000000000
mass 2 1.007940000000000
bond_coeff 1 450.000000000000000 0.957200000000000
angle_coeff 1 55.000000000000000 104.519999999999996
fix mobileEnsStr all nvt temp 300.000000000000000 300.000000000000000 100.000000000000000
fix shakeFixB1A1GwaterGrpStr all shake 0.0001 100 0 b 1 a 1
timestep 1.0
fix externalField all efield 0.0 0.0 0.1
dump dump1 all custom 1 yesEfield.dump id mol q x y z fx fy fz
dump_modify dump1 sort id
run 0
------------ getEfield.py -----------------------------
import numpy as np
def get():
#Read in data from dump files
f1 = open('noEfield.dump')
a1 = f1.readlines()[9:]
f1.close()
f2 = open('yesEfield.dump')
a2 = f2.readlines()[9:]
f2.close()
#Create arrays containing charge and force
qno = [float(i.split()[2]) for i in a1]
fzno = [float(i.split()[8]) for i in a1]
qyes = [float(i.split()[2]) for i in a2]
fzyes = [float(i.split()[8]) for i in a2]
#Compute applied electric field
eField = list()
for i in range(len(qno)):
eField.append((fzyes[i]-fzno[i])/qyes[i])
#Convert from kcal/mole-angstrom/e to V/angstrom
eField = [i*69.4786*(10**(-12)) for i in eField] #N / e
eField = [i*(1.60218*(10**19)) for i in eField] #N / C
eField = [i*(10**(-10)) for i in eField] #V/angstrom
print('max: ' + str(max(eField)))
print('min: ' + str(min(eField)))
print('mean: ' + str(np.mean(eField)))
if __name__=="__main__":
get()