gmx_ir.bsh:This script is producing NAN results when calculating the infrared spectra of TIP4P water.

[Iam-ZZG]

Dear Dr.Li, I tried the three-point water models spc, spce and the results were fine, but when calculating the infrared spectra of TIP4P water, all the data came out as NAN using gmx_ir.bsh.The version of GROMACS used is 2019.5. Below are some of the data computed. @ title "Vibrational Power Spectrum" @ xaxis label "\f{12}w\f{4} (cm\S-1\N)" @ yaxis label "a.u." @TYPE xy 0 -nan 6.67128 -nan 13.3426 -nan 20.0138 -nan 26.6851 nan 33.3564 -nan 40.0277 -nan 46.699 -nan 53.3703 -nan 60.0415 -nan 66.7128 -nan 73.3841 -nan 80.0554 -nan 86.7267 -nan 93.3979 nan 100.069 -nan 106.741 -nan 113.412 -nan 120.083 -nan 126.754 -nan 133.426 -nan 140.097 -nan 146.768 -nan 153.439 -nan 160.111 nan 166.782 -nan

[Jerkwin]

I have reproduced your findings and confirmed the problem was due to the last atom in the group you used to generate the spectra, which is the virtual site MW. Since GROMACS doesn't consider the velocities of virtual sites, the velocity autocorrelation function and spectra are also not considered. A quick fix would be to reorder the atoms in the topology, such as using OW HW1 MW HW2 instead of the original OW, HW1 HW2 MW. Another simple solution is to replace the last atom in the group used to calculate the spectra with a non-virtual site particle.

I have updated the script to use the last non-virtual site particle to store the topology, which should resolve this problem.

However, please note that I'm not certain if the approach used to calculte spectra is reliable for molecules with virtual sites. Therefore, I recommend that you carefully compare your TIP4P spectra with literature results.

[Iam-ZZG]

I am very grateful to Dr. Li for updating the script and I have tested the spectra of bulk water of TIP4P model, which is in agreement with the simulation results in the literature. But I have a question, where can I find the flexibility parameters of the water model in GROMACS? According to my search, TIP4P seems to have no flexibility model, and there is no information about the flexibility of TIP4P water model or other key parameters in TIP4P.itp in GROMACS, similarly, there is no flexibility parameter in spc.itp. I look forward to your reply, thank you very much.

[Iam-ZZG]

I noticed that there are parameters for bond energy in the ITP files, but the values of those parameters seem to be too large.

[Jerkwin]

TIP4P shares the same bonded parameters with TIP3P. You may check out the paper of TIP3P model.

[Iam-ZZG]

Thank you very much. I would like to ask if this computation method can be used to calculate the infrared spectrum of water molecules in different positions from the interface. eg.Calculating the infrared spectrum of the interface water in the figure.(I can use LAMMPS to output the total dipole moment or current vector of each water layer, and then do some post-processing).

[Jerkwin]

It should be ok, but you can not use this script directly, it needs some modifications.