GROMACS version: 2018.1
GROMACS modification: No
Greetings all,
I have been trying to run a liquid salt system with a metal cation and chloride anion using the shell MD polarization methods available in the standard GROMACS release (I have moved away from implementing Drude polarizability at this point). Unfortunately, I’m finding the documentation on shell MD to be a bit sparse.
I believe I have correctly written the topology files. The cations are not polarizable, so the anions have only the core and shell in their topology. I based the format off of information about Drude topology files and the sw.itp file recommended in an older thread. I think this file was used in the shell water study completed by GROMACS developers. I also think the topology is written correctly because GROMACS doesn’t bark at me on any steps through minimization; though, the minimization does not achieve the desired tolerance.
After energy minimization, I cannot complete an equilibration in NVT. I can clearly see that the energy and forces are large coming out of minimization. The potential energy is on the order of -1*10^6 kJ/mol and the maximum force is on the order of 10^7 kJ/mol/nm. If I check the .log file for the NVT run, it looks as though GROMACS cannot calculate even one step in the shell position energy minimization (which I have gotten it to do with a different set of polarizabilities and shell/core charges, but it never moved on to the second time step). mdrun also does not throw an error or warning when it stops running.
I think the positions of the anions and cations are okay because they come from an equilibrated simulation without polarizability. I have tried starting the shells at varying distances from their cores, somewhere between 0.001 nm and 0.3 nm away from the core and get the same results with minor changes in the energy and force. I have also attempted to tune the fcstep parameter as recommended in the manual, but I only achieve the same mdrun crash. I also have all data writing (nstxout… , nstenergy, nstlog, nstcalenergy) set to 1 so that I can see all of the output. I think my shell MD parameters are okay too (‘emtol’ = 1.0).
(1) At this point, I am not quite sure what might be causing the large energies that seem to force mdrun to stop. Is it the shell positions, fcstep, or something else?
(2) How does GROMACS calculate the spring constant? Is it based on k=(qs)^2/a from the provided core/shell charges and polarizability or is it fixed at 1000 kcal/mol/A^2 as with the Drude model and the shell water paper?
I will attach my topology, .log, and .mdp files.
Thank you for your help in advance.
polarize.top (1.8 KB)
ml.mdp (2.2 KB)
ml.log (15.1 KB)