I’m minimizing the energy in a system of ozone molecules and water solvents, and the ozone molecules are placed above the water, and separated from the water, but the energy is minimized, and the chemical bonds of the ozone molecules break. What’s going on, please? Could it be the MDP file parameters
Bonds cannot break or form in classical MD simulations or any MM process. Likely you are just seeing some kind of strained configuration that your visualization software doesn’t understand. We will need more information, including the topology, images of what you’re seeing, and .mdp settings. The output of energy minimization (maximum force and potential energy, as printed by mdrun) would also be useful.
O3-box. Gro is pre-minimization, em-vac. Gro is the result of minimization of energy. I looked at em-vac in VMD software. GRO file, where the chemical bonds between the atoms have been found to be gone. And this is my MDP file
em-vac-pme.mdp (1.3 KB)
Can you provide the topology you used for O3? It appears you do not have bonds correctly defined.
This is my topology file. Thank you.
O3.top (1.5 KB)
There are no bonds defined in the topology and you have assigned the atoms +1, +2, and +3 charges, which is physically unrealistic. Ozone should have a net zero charge. You essentially have three ions in a triangle, and you’re observing that they fly away from one another due to Coulombic repulsion.
I see. Thank you very much. But how should I adjust it, please tell me, thank you.
You will have to carefully develop suitable force field parameters by targeting QM and whatever available empirical properties you can find. Likely there are vibrational spectra for O3 for the bond and angle vibrations. Charge assignment is straightforward in a sense and is somewhat reminiscent of water. The central O will have a partial positive charge that is counterbalanced by each of the other two O atoms, which will have negative partial charges that have a magnitude half that of the central O. Determining these is somewhat challenging because there is no permanent dipole moment on the molecule, so you will have to look at higher-order multipoles (e.g. quadrupole) to compute this property and try to reproduce it with partial charge assignment.
I Got It. Thank you very much. But I still have a problem. In the previous RTP file, I defined the chemical bonds, and the amount of charge is correct (the total charge is 0) . But after converting to the top file, the bonds are gone and the charge is wrong. Why is that.
Without seeing the .rtp entry, I have no idea. Perhaps you inverted the charge and charge group columns, because pdb2gmx read the residue definition and assigned a +6 charge because it was instructed to. There are no bonds because none were likely found in the residue definition.
This is my RTP file. Is there a problem? Thank you
The partial charge and charge group columns are switched, hence why you were getting +1, +2, and +3 charges on each atom.
I see. I did write it backwards. Thank you. I would also like to ask whether the definition of my chemical bond is correct in the RTP file.
Syntactically it is correct but the two bonds should be equivalent (same force constant). I cannot assess whether the parameters are correct in the sense of giving good physical behavior.
Okay, thanks. But when I get it right, I add an aqueous solution to minimize the energy, and the ozone molecules still crack. What’s going on here？Is it because the hydrogen bonds make it look like this
Have you checked that your bond force constants are correct? They seem low to me.
Is this itIs Bond Force cosntant defined here
Yes, have you compared the force constants to other bonds in the force field?
Firstly, as Justin pointed out, there is no reason why there should be two different force constants for the two bonds. I would also expect force constants a factor 200-500 higher.
How have you determined your force constants?
I checked the information before that the two bonds of the ozone molecule is not the same, so I use the data based on two different force constants. I’m going to use the same force constant, and adjust their values to try it. Thank you
Wherever you got this information - don’t use it as a source. It’s flatly incorrect. The two bonds in ozone are equivalent. There is no formal double bond and single bond in ozone; they are equivalent via resonance.
Choosing a sensible value matters - the interplay between the vibrational frequencies of the bonds and the valence angle is important for a robust physical model.