I am trying to calculate IR spectrum of water using a flexible spc/e model. However, simple harmonic bond_style and angle_style with this model doesn’t work well. I came across a paper [1] where they have used SISM and have included some additional potential terms (Kuchitsu-Morino Intramolecular Potential) to get finer results of the vibration spectra. Unfortunately, they have not discussed how they performed these simulations.
So, first, is there any way to define integration methods (for example, SISM) other than velocity verlet algorithm for MD integrations (instead of writing the integration algorithm from scratch)?
If not, is there any way to define a custom potential term in the current verlet algorithm?
If anyone is interested, [1] J. Phys. Chem. A 2004, 108, 50, 11056–11062
If you want to implement more complex non-bonded or bonded potentials (e.g. a harmonic potential with higher order anharmonic terms added), you can modify LAMMPS to your needs. The recommended approach is to make a copy of the most similar existing potential and then rename its style and class name and make the necessary changes. There is a whole chapter in the LAMMPS manual dedicated to this: 3. Modifying & extending LAMMPS — LAMMPS documentation
In fact, most of LAMMPS was implemented by contributing authors this way. Thus, if you get it to work, we would appreciate if you contribute your adapted style(s) back to the distribution.
If you don’t want to do programming you may also be able to use a tabulated potential or use the potential styles from the LEPTON package.
Time integration in LAMMPS is done using two components: a “run style” and a time integration “fix style”. For both you can implement custom “styles”. The default setting for run_style is “verlet”. In combination with fix nve you will get the velocity verlet time stepping. There is also an alternate run style to implement r-RESPA multi-timestepping.
If you want to implement a modified Verlet, there are some options, see for example fix mvv/dpd. You can build upon that, too.
Hi Alex,
Thanks for your response. I am trying to add these anharmonic interaction in the source code and let you know how it goes.
One more doubt, is there any way to include non-bonded intramolecular potential separately? I want to keep non-bonded intermolecular interactions as they are.
One method I thought of was to define an extra bond for these interactions and keep all other potentials zero but since there is no bond among them physically, it may cause other problems.
Other method is to define a new pair potential for special bond interactions, but I couldn’t find how to define pair potential for special bonds. Can you direct me on how to include the custom potential for special bonds? Is it simply the pairwise potential (like pair_lj_cut_coul_cut.cpp)?
None of what you are asking about makes any sense to me. It sounds like you are not familiar with how LAMMPS (or any of the other MD codes I know) implement potentials.
Okay, I’ll expound a little. Let’s consider the water molecule, H1-O-H2. Now this molecule has two bonded interaction OH1 and OH2. I can use morse/ harmonic bond style to define this bond. And H1-H2 non-bonded 1-3 interaction. To define the potential due to this interaction I can use special bonds. This special bond potential is a custom potential I want to define. That’s my question, how to define this potential for a special bond?
You are still not making sense.
If the water model is a rigid model, the H1-H2 interaction is ignored and thus irrelevant.
If the model is a flexible model, the H1-H2 interaction will be covered by the H1-O-H2 angle potential and any non-bonded contribution should be excluded to avoid unphysical double counting of that interaction.
That’s the point. I am not trying to simulate the rigid water model. I want to get the IR spectra as I said earlier. Which depends on the intramolecular vibrations. My water model is flexible SPC/E. and the interaction between H1-H2 depends on the bond lengths of OH1, OH2 and the angle H1-O-H2. This is in addition to the H1-O-H2 angle potential.
As I already stated, this must be implemented as part of the angle potential. See for example the Urey-Bradley terms in the CHARMM force field angle interactions: angle_style charmm command — LAMMPS documentation
Ohh right. I can include it in the angle potential itself. For some reason I was stuck in the bond potential. Thanks, this will help. (Y)
The term you describe is neither a bond potential nor a non-bonded pairwise potential, but a manybody interaction term that requires knowing all three interacting particles in the same molecule, which are two bonds sharing one atom. Covering this with an angle potential is the only logical solution.
1 Like