According to what you said, I finally simulated successfully. But I still have a question:
In fact, I want to use SPC/E water model instead of SPC water model. So I try to only use the spce.lt file come with Moltemplate as the parameters of force field of water and use oplsaa.lt file as parameters of force field of others. However, it failed with the same error. So I guess is it because the SPC/E water model can’t be applied to OPLS-AA force field?
There is indeed something wrong when I use the SPC/E model of water, which I suspect is related to the topology that Moltemplate creates when mixing a model with the bonds defined as a list (the polymer) and one with the bonds explicitly assigned (the water).
In principle, I wouldn’t expect a big difference as the only difference is with the charges. Give me a few days to find the error in the topology, then I’ll try to fix it.
FWIW, OPLS-AA has been parameterized for use with TIP3P or (by preference) TIP4P. CiteSeerX
Absolutely. Plus, mixing force fields with different parametrization strategies is never a good idea.
There is probably a bug in Moltemplate, as the lost atoms point to something fishy with the topology.
EDIT:
The initial geometry was to blame.
The problem is bad contact between water molecules no. 12 and 270. Somehow, the minimisation fixes the geometry with SPC charges but not with SPC/E charges, which are slightly bigger.
The solution is either to do a quick annealing of the system and then start a new simulation with SPC/E water, or to use a different minimisation algorithm. Here is a minimum example that work with your initial geometry and SPC/E water:
spce.lt
# This is a version of the SPC/E water model
# that goes with the OPLS-AA force field.
SPCE {
write("Data Atoms") {
$atom:O $mol:. @atom:O -0.8476 0.0000000 0.00000 0.000000
$atom:H1 $mol:. @atom:H 0.4238 0.8164904 0.00000 0.5773590
$atom:H2 $mol:. @atom:H 0.4238 -0.8164904 0.00000 0.5773590
}
write_once("Data Masses") {
@atom:O 15.9994
@atom:H 1.008
}
write("Data Bonds") {
$bond:OH1 @bond:OH $atom:O $atom:H1
$bond:OH2 @bond:OH $atom:O $atom:H2
}
write("Data Angles") {
$angle:HOH @angle:HOH $atom:H1 $atom:O $atom:H2
}
write_once("In Settings") {
bond_coeff @bond:OH 600.0 1.0
angle_coeff @angle:HOH 75.0 109.47
pair_coeff @atom:O @atom:O 0.1553 3.166
pair_coeff @atom:H @atom:H 0.0 0.0
}
}
sample04.lt
# Use the OPLS-AA force field for all species except water, which is SPC/E.
import "oplsaa.lt"
import spce.lt
import "PolyNIPAM.lt"
# Define the SPC water and ions as in the OPLS-AA
Ca inherits OPLSAA {
write("Data Atoms"){
$atom:a1 $mol:. @atom:354 0.0 0.00000 0.00000 0.000000
}
}
Cl inherits OPLSAA {
write("Data Atoms"){
$atom:a1 $mol:. @atom:344 0.0 0.00000 0.00000 0.000000
}
}
# Create the system.
wat=new SPCE[500]
pol=new PolyNIPAM[1]
cat=new Ca[1]
ani=new Cl[2]
# Periodic boundary conditions:
write_once("Data Boundary"){
0 26 xlo xhi
0 26 ylo yhi
0 26 zlo zhi
}
# Define the input variables.
write_once("In Init"){
# Input variables.
variable ts equal 1 # timestep
variable temp equal 298.15 # equilibrium temperature
variable p equal 1. # equilibrium pressure
# PBC (set them before the creation of the box).
boundary p p p
neighbor 3 bin
}
# Let's go, Pickacku!
write_once("In Run"){
compute pe1 all pe/atom pair
dump TRJ all custom 100 sample04.dump id xu yu zu c_pe1
timestep \$\{ts\}
thermo 1000
thermo_style custom step etotal evdwl ecoul elong ebond eangle edihed eimp pe ke temp press atoms vol density
thermo_modify flush yes
min_style fire/old
minimize .001 .001 1000 100000
write_data sample04.min
group spce type @atom:O @atom:H
fix fRattleSPCE spce rattle 0.0001 10 100 b @bond:OH a @angle:HOH
# Reset the velocity after the minimisation, to get rid of excess kinetic energy.
velocity all create \$\{temp\} 12345
fix 1 all nvt temp \$\{temp\} \$\{temp\} \$(100*dt) drag 2
run 2000
unfix 1
}
The topology generated by Moltemplate was indeed correct.
Thanks for your advice and I’ll try to use a different minimization algorithm.
And now I have one more question: If I can’t assign types for every atoms in molecules because there isn’t such type in oplsaa.lt file, I could only choose similar types but the result of doing this is that the charge of the entire molecule is not what I expected. Is it because that the type I assigned is incorrect？I wonder whether I could adjust the charge appropriately on my own.
In this case, you must parametrize the molecule using the same approach as the original OPLS-AA force field. In other words, the Lennard-Jones parameters and partial charges must be consistent with OPLS-AA, otherwise, you will likely learn very little from the simulation. We assume we will use existing bond, angle, dihedral, and improper terms in the OPLS-AA force field, but we specify new charges and pair_coeff terms.
It shouldn’t be difficult to pass those parameters to Moltemplate, though. Just create a separate LT file for your molecule that contains the following information:
- The macro defining the new molecule will inherit the bonded parameters (e.g.
NEWMOL inherits OPLSAA {}). - In the
write('Data Atoms') {}block, use new atom types (e.g.@atom:myC) and assign the charges directly in this block. - Specify the connectivity only with the
write('Data Bond List') {}block. - Specify the new LJ coefficients:
write_once("In Settings") { pair_coeff @atom:myC @atom:myC 0.066 3.5 ... } - Use Moltemplate’s
replacecommand to rename your new atom types, adding explicit bonded interaction parameters to the name: e.g.
At this stage, you decide which bonded coefficient you want for your new atom types. Work out how to rename the new atom types by looking for similar ones in the filereplace{ @atom:myC @atom:001_b013_a013_d013_i013 } # alkyl C in CH3 replace{ @atom:myH @atom:002_b046_a046_d046_i046 } # alkyl H in CH3oplsaa.lt.
Alternatively, you will need to explicitly specify all the bonded and non-bonded parameters for your molecule. In my experience, bonded interactions do not differ too much in type I force fields. If you want to increase the accuracy of your model, you can specify custom dihedral potentials, as they are more sensitive to the molecular topology and hence more difficult to generalise.
Good luck, and let us know how your project goes.
If I still use SPC/E water, even if I try to use minimize style ‘fire’ (minimize style ‘fire/old’ has been removed from LAMMPS after the 22 December 2022 version.), it’ll generate a new error.
As I have explained, the error is due to overlapping molecules in your input geometry that the minimiser cannot disentangle when using SPC/E charges. Just use the SPC model as in my first working example, save a simple DUMP with XYZ coordinates, and restart the simulation from there, switching to SPC/E water.
Is it so difficult to understand why your simulation keeps crashing?
Thanks, now I understand.
Is there a software you know of that can translate a molecule in .xyz or .pdb to a LAMMPS topology file in the OPLS-aa? Or is it better to just brute force it ourselves?
If you haven’t noticed, this thread is about describing this very task using Moltemplate.
A more didactic How-To based on the example described here is now present in the latest version of the ManualTM, section 8.6.4 Moltemplate Tutorial.
1 Like
Do you think, given this post, moltemplate can properly handle that?
Of course it can. However, the version of OPLS-AA shipped with Moltemplate is a bit outdated, as highlighted in this Question regarding Napthalene C (atom type 92) in OPLSAA · Issue #93 · jewettaij/moltemplate · GitHub.
You can still use it to define a new molecule and, as long as the atom types are assigned correctly, you will have a compliant OPLS molecular topology. Charge neutrality is a good indicator.
1 Like