Velocity Autocorrelation Fucntion for SPC/E water model

LAMMPS LAMMPS Beginners
1

Hello.

I wanted to calculate VACF for the SPC/E water model where in the input scripts of LAMMPS (attached below), I read the data of the coordinates of all the atoms (9576 atoms total for 3192 molecules) from the GW.dat file and got the values of VACF using “compute” command. However, I got a similar trend of VACF but did not get the accurate curve. I have just started to learn MD and am still a novice, so can anyone tell me how to improve my result even more by looking at the script or give any advise regarding that? Thank you.

###################################
# Simulation settings
###################################
units           real
atom_style      full
boundary	p p p

####################################
# Reading data
####################################
read_data       GW.dat
group		water type 1 2
group		oxygen type 2
group		hydrogen type 1

#####################################
# Force interactions
#####################################
pair_style  	lj/cut/coul/long 10
bond_style      harmonic
angle_style     harmonic
kspace_style    pppm 1.0e-6

pair_coeff 1 1 0.0 0.4000 #H-H
pair_coeff 2 2 0.15535 3.1660 #O-O
pair_coeff 1 2 0.0 1.7753 #H-O

bond_coeff	1 1000 1 # H-O
angle_coeff	1 1000 109.47 # H-O-H

#################################################
# Variables
#################################################
variable 	Nrun equal 100000
variable 	Nf equal ${Nrun}/100
variable 	Ne equal 10
variable 	Nr equal ${Nf}/${Ne}
variable 	Ndump equal ${Nrun}/2
variable 	Nr_rdf equal 0.5*${Nrun}/${Ne}

variable 	watMoleMass equal 18.0153 # /(g/mol)
variable 	nAvog equal 6.0221415e23 # Avogadro's number
variable 	watMoleculeMass equal (${watMoleMass}/${nAvog}) # /(g/molecule)
variable 	A3_in_cm3 equal 1e-24 # Angstrom^3 in cm^3
variable 	nAtoms equal atoms
variable 	nMolecules equal v_nAtoms/3

variable 	Text equal 298.0
variable 	Pext equal 1.0

#################################################
# Initial conditions and settings
#################################################
velocity 	water create ${Text} 1234546 dist gaussian
neighbor 	2.0 bin
neigh_modify 	every 1 delay 0 check yes

##############
# Minimization
##############
minimize        1.0e-5 1.0e-7 1000 10000

##################################################
# Computations
##################################################
reset_timestep	0
timestep	1

#----------- NPT Equilibration -------------------
fix 		SHAKE water shake 1e-6 200 0 b 1 a 1
fix             NPT water npt temp ${Text} ${Text} 100 iso ${Pext} ${Pext} 1500
#fix 		removeMomentum water momentum 1 linear 1 1 1

compute 	T water temp
fix 		TempAve water ave/time ${Ne} ${Nr} ${Nf} c_T 
variable 	P equal press
fix 		PressAve water ave/time ${Ne} ${Nr} ${Nf} v_P
compute 	PE all pe pair kspace
variable 	PE_Mol equal c_PE/v_nMolecules
fix 		PEAve_Mol all ave/time ${Ne} ${Nr} ${Nf} v_PE_Mol
variable 	Dens equal v_nMolecules*${watMoleculeMass}/(vol*${A3_in_cm3})
fix 		DensAve water ave/time ${Ne} ${Nr} ${Nf} v_Dens file wat.dens

thermo_style 	custom step temp f_TempAve press f_PressAve f_PEAve_Mol f_DensAve 
thermo_modify 	flush yes
thermo		1000


#------------- NVE equilibration ------------------
unfix           NPT
#unfix 		removeMomentum
fix             NVE water nve

compute vac all vacf
fix myvcf all ave/time 2 1 20 c_vac[1] c_vac[2] c_vac[3] c_vac[4] file vcf2.dat


run		1000
2

It might help if you tell us how your results compared with what you expected. 1000 steps is pretty short, and you might need longer runs for more accurate results. Your input script is long enough that I can’t obviously notice anything wrong. (This might mean there is nothing wrong, or that there is something wrong but it’s not obvious.)

3

Thank you for your reply.

This is the figure that I got (Normalized VACF vs Timestep):
VCFN

But I expected this type of figure:

vcfp
vcfp700×525 24.4 KB

Is there be problem with frequency or NVE run?

4

I could just say “your run isn’t long enough”. But I am just a stranger on the Internet and you will gain much more from my advice if you ask yourself: “How would I determine that what srtee said was correct or not?”

5

After running 1000 times in NPT and 100,000 times in NVE and computing vacf under NVE, I got this:

image
image700×525 57.5 KB

1 Like
6

I would also suggest that after you minimize the initial structure, you should probably do an equilibration run before grabbing data in a production run. But, as @srtee says, you need to understand the reasons behind what you do, because you can’t write in your paper “I did this because someone on the internet said to”.

So you need to find answers for “should I use a minimized structure directly in a production run?” and “what is a reasonable time for accumulating data?” and “will one reference point (the beginning of the run) be enough statistically for my system?”.