units real
boundary p p p
atom_style full
read_data data.glycine_water
# define groups [http://lammps.sandia.gov/doc/group.html]
group glycine type 1 2 3 4 6 7 #defines the group glycine.
group water type 5 8 #defines the group water.
# interaction styles
# pair_styles documentation page: http://lammps.sandia.gov/doc/pair_lj.html
pair_style lj/cut/coul/long 12.0 12.0
kspace_style pppm 1.0e-4
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_modify mix geometric tail no
# tail means long range correction factors added to energy/pressure. mix command defined below:
# mix geometric
# epsilon_ij = sqrt(epsilon_i * epsilon_j)
# sigma_ij = sqrt(sigma_i * sigma_j)
#
# mix arithmetic
# epsilon_ij = sqrt(epsilon_i * epsilon_j)
# sigma_ij = (sigma_i + sigma_j) / 2
# OPLS considers 1-4 interactions with 50%.
special_bonds lj/coul 0.0 0.0 0.5
# force field parameters
# missing nonbonded parameters are inferred from mixing.
#field, atom type, atom type, eps(kcal/mol), sigma(angstroms) [http://lammps.sandia.gov/doc/pair_coeff.html]
# 1 C1
# 2 C2
# 3 H1
# 4 H2
# 5 HW
# 6 N
# 7 O
# 8 OW
pair_coeff 1 1 0.1050 3.7500 # C1
pair_coeff 2 2 0.1050 3.5000 # C2
pair_coeff 3 3 0.0000 0.0000 # H1
pair_coeff 4 4 0.0300 2.5000 # H2
pair_coeff 5 5 0.0000 0.0000 # HW
pair_coeff 6 6 0.1700 3.2500 # N
pair_coeff 7 7 0.2100 2.9600 # O
pair_coeff 8 8 0.1554 3.1654 # OW
#field, bond type, kb, r0 [http://lammps.sandia.gov/doc/bond_coeff.html]
# 1 C1-C2
# 2 C1-O
# 3 C2-H2
# 4 C2-N
# 5 H1-N
# 6 HW-OW
bond_coeff 1 317.00 1.5220 # C1-C2
bond_coeff 2 656.00 1.2500 # C1-O
bond_coeff 3 340.00 1.0900 # C2-H2
bond_coeff 4 367.00 1.4700 # C2-N
bond_coeff 5 434.00 1.0100 # H1-N
bond_coeff 6 1059.162 1.000 # HW-OW
#field, angle type, kb, theta0 [http://lammps.sandia.gov/doc/bond_coeff.html]
# 1 C1-C2-H2
# 2 C1-C2-N
# 3 C2-C1-O
# 4 C2-N-H1
# 5 H1-N-H1
# 6 H2-C2-H2
# 7 H2-C2-N
# 8 HW-OW-HW
# 9 O-C1-O
angle_coeff 1 35.0 109.5 # C1-C2-H2
angle_coeff 2 80.0 111.2 # C1-C2-N
angle_coeff 3 70.00 117.00 # C2-C1-O
angle_coeff 4 35.0 109.5 # C2-N-H1
angle_coeff 5 35.0 109.5 # H1-N-H1
angle_coeff 6 33.0 107.80 # H2-C2-H2
angle_coeff 7 35.0 109.50 # H2-C2-N
angle_coeff 8 75.9 113.24 # HW-OW-HW
angle_coeff 9 80.0 126.0 # O-C1-O
#field, dihedral type, K1, K2, K3, K4 [http://lammps.sandia.gov/doc/dihedral_coeff.html]
# 1 C1-C2-N-H1
# 2 H2-C2-N-H1
# 3 O-C1-C2-H2
# 4 O-C1-C2-N
dihedral_coeff 1 0.000 0.000 0.347 0.000 # C1-C2-N-H1
dihedral_coeff 2 0.000 0.000 0.300 0.000 # H2-C2-N-H1
dihedral_coeff 3 0.000 0.000 0.000 0.000 # O-C1-C2-H2
dihedral_coeff 4 0.000 0.546 0.000 0.000 # O-C1-C2-N
#field, improper type, K1, X0 [http://lammps.sandia.gov/doc/improper_coeff.html]
improper_coeff 1 7.000 180.0 #C2-O-C1-O
run_style verlet #The verlet style is a standard velocity-Verlet integrator.
# velocity command documentation: http://lammps.sandia.gov/doc/velocity.html
velocity all create 300.0 239847 dist gaussian mom no rot no
# the rundown:
# "velocity all" prompts lammps that it will generate random velocities for all molecules in the simulation
# "create 300 239847" is the actual command to create random velocities. 300 is the temp (Kelvin) and 239847 is the seed
# "dist gaussian" makes the velocity distribution gaussian as opposed to uniform
# "mom no" prevents lammps from zeroing all linear momenta
# "rot no" does the same thing for rotational momenta
# fix command documentation: http://lammps.sandia.gov/doc/fix.html
fix 1 all nvt temp 300.0 300.0 100.0
# the rundown:
# "fix 1 all" --2 identifies this instance of the command, all identifies all particles affected
# "nvt temp" makes a thermostat. See documentation for the particulars.
# "300 300 100.0" means Tstart, Tstop, damping parameter
# Optional, as far as I can tell: change the timestep to a smaller value. LAMMPS real units are femtoseconds.
timestep 0.25 #resets timestep to 0.25 fs
dump 98 all xyz 100 TTA-start.xyz
# run just a few steps to get rid of all the symmetry
run 100
# minimize command documentation: http://lammps.sandia.gov/doc/minimize.html
minimize 1.0e-5 1.0e-6 10000 100000
# the rundown:
# "1.0e-5" stopping tolerance for energy (unitless)
# "1.0e-6" stopping tolerance for force (force units)
# "10000" max iterations of minimizer
# "100000" max number of force/energy evaluations
# Write a restart file after minimization, which can be picked up and used as a checkpoint later if necessary
write_restart glycine_water-min.restart
reset_timestep 0 #does not reset the timestep between force calculations, but changes the counter to 0. Do this after minimization
# neigh_modify command documentation: http://lammps.sandia.gov/doc/neigh_modify.html
neigh_modify every 5 delay 0 check no
# the rundown:
# "every 5" updates the verlet neighbor list every 1 timesteps
# "delay 0" only build every timestep. What is "build"?
# "check no" build even if no atom has moved half the skin distance or more
neighbor 3.0 bin
# command: neighbor skin style
# skin = skin thickness in angstroms for real units (default is 2)
# style = calculation method. "bin" scales linearly, "nsq" scales quadratically, so always use "bin"
thermo 20 #compute and print thermodynamic information every 20 timesteps.
thermo_style multi #prints a multiple-line summary of information. Can also use "one" or "custom ___" where ___ is other arguments
velocity all scale 300 #scales all particle velocities to accomplish a temperature of 300 K
# dump command documentation: http://lammps.sandia.gov/doc/dump.html
dump 1 all dcd 20 glycine-water-eq.dcd
# the rundown:
# "dump 1 all" -- 1 is the user assigned name for the dump, just like "fix" identifiers. "all" means all particles
# "dcd" -- the type of file to save it as
# "20" means dump every 20 timesteps
# "glycine-water-eq.dcd" is the name of the file to save it as
dump 2 all xyz 20 glycine_water_min.xyz
# dump_modify command documentation: http://lammps.sandia.gov/doc/dump_modify.html
dump_modify 1 unwrap no
# the rundown:
# "1" --another user defined identifier for the command
# "unwrap no" --a command about how to store the data, only for .dcd or .xtc files. It means the atom positions will NOT appear ...
# ... as if there were no periodic boundary conditions in the saved file.
compute rdf1 all rdf 200 7 1 7 3 7 4 7 5 6 1 3 4 4 5 4 4 5 5 6 8
# the rundown:
# "rdf1" is the name of the variable in which data is saved
# "all" means all atoms
# "200" means it divides into 200 bins
# remaining numbers are the pairs between which the rdf is computed.
# 7 1 is O - C1 (O in glycine to C1(carboxylate carbon) in glycine)
# 7 3 is O - H1 (O in glycine to H1(on nitrogen) in glycine)
# 7 4 is O - H2 (O in glycine to H2(on carbon) in glycine)
# 7 5 is O - HW (O in glycine to HW in water)
# 6 1 is N - C1 (N in glycine to C1(carboxylate carbon) in glycine)
# 3 4 is H1 - OW (H1(on nitrogen) in glycine to OW in water)
# 4 5 is OW - HW (water oxygen to water hydrogen)
# 4 4 is OW - OW (water oxygen to water oxygen)
# 5 5 is HW - HW (water hydrogen to water hydrogen)
# 6 8 is N - OW (N in glycine to O in water)
compute wat_diff water msd
compute gly_diff glycine msd
fix 2 all ave/time 10 800 12000 c_rdf1 file rdf1.rdf mode vector
# fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ...
# ID, group-ID are documented in fix command
# ave/time = style name of this fix command
# Nevery = use input values every this many timesteps
# Nrepeat = # of times to use input values for calculating averages
# Nfreq = calculate averages every this many timesteps
# one or more input values can be listed
# value = c_ID, c_ID[N], f_ID, f_ID[N], v_name
# For example, if Nevery=2, Nrepeat=6, and Nfreq=100, then values on timesteps 90,92,94,96,98,100
# will be used to compute the final average on timestep 100.
fix 3 all ave/time 100 1 100 c_wat_diff file diff_wat.txt mode vector
fix 4 all ave/time 100 1 100 c_gly_diff file diff_gly.txt mode vector
#
run 24000
#save a restart file of this as well
write_restart gly_wat_eq.restart
