Skip to content

Commit

Permalink
Fixed all the water models (SPC, SPC/E, TIP3P, TIP5P). They should al…
Browse files Browse the repository at this point in the history
…l be working again with OPLSAA (2023 version). I alwo created examples for each of these water models and tested them. (They are all working, but there is a problem with the NVT protocol in the TIP5P example. Since NPT is working for this example, fixing NVT is a lower priority. I'll fix this later when I have time.)
  • Loading branch information
jewettaij committed Dec 7, 2024
1 parent 9dea361 commit f027da1
Show file tree
Hide file tree
Showing 40 changed files with 3,212 additions and 115 deletions.
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
import "ethylene.lt" # <- defines the "Ethylene" molecule type.
import "benzene.lt" # <- defines the "Benzene" molecule type.


# Periodic boundary conditions:
write_once("Data Boundary") {
0.0 48.00 xlo xhi
Expand Down

This file was deleted.

Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import "spce_oplsaa.lt" # <- defines the "SPCE" (water) molecule type.
import "methane.lt" # <- defines the "Methane" molecule type (uses OPLSAA)
import "spce_oplsaa.lt" # <- defines the "SPCE" (water) molecule type.

# Periodic boundary conditions:
write_once("Data Boundary") {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ write_data "system_after_min.data"


timestep 1.0
dump 1 all custom 5000 traj_npt.lammpstrj id mol type x y z ix iy iz
dump 1 all custom 2000 traj_npt.lammpstrj id mol type x y z ix iy iz


# Turn on the barostat and thermostat
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
The purpose of this example is to test the density of water
constructed using the OPLSAA force-field. (I think this is SPC water, not SPCE)

I just wanted some kind of sanity check to make sure we are converting
the OPLSAA parameters into moltemplate/LAMMPS format correctly.

The "TEST_density_estimate.txt" contains the results of that test.

-------- Instructions: ---------

More detailed instructions on how to build LAMMPS input files and
run a short simulation are provided in other README files.

step 1)
README_setup.sh

step 2)
README_run.sh


### Customizing atomic charges

In this example, atomic charge for OPLSAA atoms is determined by @atom type
(...according to a lookup table located at the beginning of the
"oplsaa.lt" file).
(Note: Any atomic charges listed in the "Data Atoms" section will be ignored.)
**These charges can be overridden.**
See the "README.md" file located in the parent directory
for instructions explaining how to customize atomic charge.
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@

# Note: By default, the system.data and system.in.settings files contain
# extra information for atoms defined in OPLSAA which you are not using
# in this simulation.
# This is harmless, but if you to delete this information from your
# system.in.settings and system.in.data files, run this script:

cleanup_moltemplate.sh

# (Note: Removing unecessary atom types will make it easier to visualize the
# simulation in VMD.)
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
# --- Running LAMMPS ---
# The 2 files "run.in.npt", and "run.in.nvt" are LAMMPS
# input scripts which link to the input scripts and data files
# you hopefully have created earlier with moltemplate.sh:
# system.in.init, system.in.settings, system.data
# If not, carry out the instructions in "README_setup.sh".
#
# -- Instructions: --
# If "lmp_mpi" is the name of the command you use to invoke lammps,
# then you would run lammps on these files this way:


lmp_mpi -i run.in.npt # minimization and simulation at constant pressure
lmp_mpi -i run.in.nvt # simulation at constant volume


# If you have compiled the MPI version of lammps, you can run lammps in parallel
#mpirun -np 4 lmp_mpi -i run.in.npt
#mpirun -np 4 lmp_mpi -i run.in.nvt
# (assuming you have 4 processors available)
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@

# Create LAMMPS input files this way:
cd moltemplate_files

# run moltemplate

# This was the original (simple) way to run moltemplate:
# moltemplate.sh system.lt <-- COMMENTING OUT
# Instead, this is the recommended way to run moltemplate with OPLSAA:

moltemplate.sh system.lt -report-duplicates bytype __

# (The optional "-report-duplicates bytype __" arguments check to make
# sure that there was no ambiguity in the dihedrals that were generated.
# This is an issue with OPLSAA. If there was, then moltemplate will create
# a file named "warning_duplicate_dihedrals.txt".)
#
# (Note: You can also check for missing angle,dihedral params this way:)
# moltemplate.sh -checkff system.lt


# Moltemplate generates various files with names ending in *.in* and *.data.
# Move them to the directory where you plan to run LAMMPS (in this case "../")
mv -f system.data system.in* ../

# Optional:
# The "./output_ttree/" directory is full of temporary files generated by
# moltemplate. They can be useful for debugging, but are usually thrown away.
rm -rf output_ttree/

# Optional:
rm -f run.in.EXAMPLE # not needed. We have several run.in... files already.

# Optional:
# If any warnings were generated, move them to the parent folder
# (so they get noticed).
mv -f warning*.txt ../ 2> /dev/null

cd ../




# Optional:
# Note: The system.data and system.in.settings files contain extra information
# for atoms defined in OPLSAA which you are not using in this simulation.
# This is harmless, but if you to delete this information from your
# system.in.settings and system.in.data files, run this script:
#
# cleanup_moltemplate.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@

------- To view a lammps trajectory in VMD --------


1) Build a PSF file for use in viewing with VMD.

This step works with VMD 1.9 and topotools 1.2.
(Older versions, like VMD 1.8.6, don't support this.)


a) Start VMD
b) Menu Extensions->Tk Console
c) Enter:

(I assume that the the DATA file is called "system.data")

topo readlammpsdata system.data full
animate write psf system.psf

2)

Later, to Load a trajectory in VMD:

Start VMD
Select menu: File->New Molecule
-Browse to select the PSF file you created above, and load it.
(Don't close the window yet.)
-Browse to select the trajectory file.
If necessary, for "file type" select: "LAMMPS Trajectory"
Load it.

---- A note on trajectory format: -----
If the trajectory is a DUMP file, then make sure the it contains the
information you need for pbctools (see below. I've been using this
command in my LAMMPS scripts to create the trajectories:

dump 1 all custom 5000 DUMP_FILE.lammpstrj id mol type x y z ix iy iz

It's a good idea to use an atom_style which supports molecule-ID numbers
so that you can assign a molecule-ID number to each atom. (I think this
is needed to wrap atom coordinates without breaking molecules in half.)

Of course, you don't have to save your trajectories in DUMP format,
(other formats like DCD work fine) I just mention dump files
because these are the files I'm familiar with.

3) ----- Wrap the coordinates to the unit cell
(without cutting the molecules in half)

a) Start VMD
b) Load the trajectory in VMD (see above)
c) Menu Extensions->Tk Console
d) Try entering these commands:

pbc wrap -compound res -all
pbc box

----- Optional ----
Sometimes the solvent or membrane obscures the view of the solute.
It can help to shift the location of the periodic boundary box
To shift the box in the y direction (for example) do this:

pbc wrap -compound res -all -shiftcenterrel {0.0 0.15 0.0}
pbc box -shiftcenterrel {0.0 0.15 0.0}

Distances are measured in units of box-length fractions, not Angstroms.

Alternately if you have a solute whose atoms are all of type 1,
then you can also try this to center the box around it:

pbc wrap -sel type=1 -all -centersel type=2 -center com

4)
You should check if your periodic boundary conditions are too small.
To do that:
select Graphics->Representations menu option
click on the "Periodic" tab, and
click on the "+x", "-x", "+y", "-y", "+z", "-z" checkboxes.

5) Optional: If you like, change the atom types in the PSF file so
that VMD recognizes the atom types, use something like:

sed -e 's/ 1 1 / C C /g' < system.psf > temp1.psf
sed -e 's/ 2 2 / H H /g' < temp1.psf > temp2.psf
sed -e 's/ 3 3 / P P /g' < temp2.psf > system.psf

(If you do this, it might effect step 2 above.)
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
This system contains 1728 water molecules
(This is SPC water I think.)

Then I ran a short simulation for 170000 timesteps at 300Kelvin and 1 atm.
(that's when it crashed. I'll worry about why later...)

Anyway, the average volume was 52149.8 (in Angstroms^3)
(for the last 80000 timesteps, after it had equilibrated)

Given that the mass of water is 18.0154 grams per mole, I'm getting
this value for the density:

density = (1728*18.0154/6.02214129e23) / (52149.8*1e-30*1e6)
= 0.991 (in grams per mL)

I'm only looking for gross errors in the OPLSAA force-field.
So I'm satisfied with a 1% error.
But I realize this is not a particularly rigorous test.

Andrew 2014-5-21
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
import "spc_oplsaa.lt" # <- defines the "SPC" (water) molecule type

# NOTE: The "spc_oplsaa.lt" file is located in the moltemplate/force_fields/
# folder (distributed with moltemplate).

# Periodic boundary conditions:
write_once("Data Boundary") {
0.0 41.40 xlo xhi
0.0 41.40 ylo yhi
0.0 41.40 zlo zhi
}

# The next command generates a (rather dense) cubic lattice with
# spacing 3.45 Angstroms. (The pressure must be equilibrated later.)

waters = new SPC [12].move(0.00, 0.00, 3.45)
[12].move(0.00, 3.45, 0.00)
[12].move(3.45, 0.00, 0.00)

Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
# PREREQUISITES:
#
# You must use moltemplate.sh to create 3 files:
# system.data system.in.init system.in.settings
# (Follow the instructions in README_setup.sh,
# or run the file as a script using ./README_setup.sh)
#
# ------------------------------- Initialization Section --------------------

include "system.in.init" # specify the style of force field used

# ------------------------------- Atom Definition Section -------------------

read_data "system.data"

# ------------------------------- Settings Section --------------------------

include "system.in.settings" # load the force field parameters
include "system.in.charges" # load the charge of each atom

# ------------------------------- Run Section -------------------------------


# -- minimization protocol --

thermo 50
minimize 1.0e-4 1.0e-6 100000 400000

unfix fShakeSPC # Disable SHAKE during minimization and pressure equilibr

# Optional: write the coordinates after minimization
write_data "system_after_min.data"


# -- simulation protocol --


timestep 1.0
dump 1 all custom 1000 traj_npt.lammpstrj id mol type x y z ix iy iz
fix fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0
thermo 100

run 100000

write_data "system_after_npt.data"
Loading

0 comments on commit f027da1

Please sign in to comment.