JustPaste.it

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