Some benchmark tests are required. Often this is done when you write a proposal.

I'm currently running solute tempering on 2 nodes (or 48 cores) per replica. The only difference I can see between the setups - also related to the error message you attached - is the number of OMP threads per MPI process: in my case it's 1 thread per process, in yours it's 2. You can specify this ratio manually using the -ntomp option of gmx_mpi mdrun, try and see if that fixes the issue.

Hi, Miłosz Wieczór ,

Could you let me know what command are you using for submission? The command I used was:

mpirun -np 60 gmx_mpi mdrun -ntomp 1 -deffnm topol -plumed ../plumed.dat -multi 20 -replex 200 -hrex -cpi topol

and the output I got is attached. It says the system is not equilibrated enough, but if I used the exactly same input files and just change the submission command to 2 cores per replica, the simulation works fine. Could you let me know if you have any suggestion? I greatly appreciate your help!

Regards,

Peiyin

The thing is that when the simulation goes awry, LINCS and inter-process communication are the first to report an error. If you have neither (as would happen with 2 cores, where 1 is assigned to PME and another one to real-space interactions), the job might survive the disturbance without errors, but that would likely be visible in the trajectory. Hence more cores = less tolerance for rapid jumps.

If you monitor any reaction coordinates in Plumed, you might temporarily use very low (=1) STRIDE and FLUSH to see if anything suspicious happens to the structure. Otherwise you could use a low value of nstxout-compressed in .mdp to see where the problem might be.

You can also try out my script to prepare Gromacs topologies for solute tempering runs to see if there might be the problem: https://gitlab.com/KomBioMol/solute_tempering

Hi, Miłosz Wieczór ,

Thank you very much for your reply. It is very instructive.

May I ask what is the time step you use for your simulation? I have reduced the time step to 0.5 fs from 2 fs but still see issue sometimes.

I don't monitor any collective variables through Plumed. My plumed file is empty.

I will try reducing the frequency of output file to 1 and see if I can observe anything that is wrong.

Also, I didn't explain my procedure very clearly in the previous posts. The .tpr files I generated for REST2 used .gro files after the NPT equilibration so it should be pretty stable. I can run regular MD production simulation using this .gro file with 2 fs and everything can be run smoothly. However, if I use REST2 (the -hrex flag), I can not run it with 2 fs and the error will say that the system is not well equilibrated enough, so that's why I'm reducing the time step to 0.5 fs. I don't know if it is because I'm using too many replicas

(20 in this case) and thus increasing the possibility of error occurrence. Do you have any suggestion?

Thank you again!

I saw that you chose 1 fs for the run that failed, which made me suspicious about LINCS - if you're using it, you should at least get LINCS warnings in the logfiles from the 2-thread runs. You might check if that is the case. By the way, what is the highest temperature that you use for your replicas?

I also noticed the instability doesn't occur at the very beginning of your simulation, prompting me to question the topology/parameter files. I use the typical 2 fs timestep with 12 replicas and had no issues over more than 30 ns now - the number of replicas is clearly not an issue. My wild guess here would be an incorrectly rescaled force field parameter as a culprit, which is why I suggested replicating the issue with independently generated topologies.

Hi Miłosz Wieczór ,

Yes, I am using LINCS. I have checked the .log files and only saw the following thing that is related to LINCS, which doesn't look like a warning to me:

P-LINCS: A Parallel Linear Constraint Solver for molecular simulation

Estimated maximum distance required for P-LINCS: 0.218 nm

The temperature range I use is 300-500K.

May I ask how do you generate the .top files and the .tpr files? I follow the instructions on the Plumed tutorial (https://www.plumed.org/doc-v2.4/user-doc/html/hrex.html) and use their commands. Also, I have posted another question that's related to my system:

https://www.researchgate.net/post/REST2_Infinite_potential_fatal_error_occurs_when_REST2_tries_to_swap_why

In that question I have attached my .mdp file and the error output file for a new run I just did using a .gro file from a 360-ns long MD simulation so the system should be quite stable in this case. Since the error occurs at step 201, which is the exchange frequency (200) I set for REST2, I highly suspect the reason for the failure is due to the Hamiltonian exchange swap. However, I don't know why the simulation will end up with infinite potential when 2 configurations swap. I greatly appreciate your help as I have been stuck with this issue for a long while. Thank you so much!

Hi, Miłosz Wieczór ,

I was able to generate .top files using your Python scripts but couldn't generate the .tpr files because the format is different (incorrect number of parameters, please see the attached screenshot). However, I compared the content of the .top files from the ones I generated with your scripts and the ones generated using Plumed's "partial_tempering" command and they are basically the same, with just different notations when specifying the "heated atoms". The charges scaled are the same using two ways, so I think my .tpr files should be fine in terms of the scaling potential...

Ok, I see - indeed, my script does not yet support the Ryckaert-Belleman dihedrals used in OPLS. But since you used plumed's scripts, and the values are similar, there should be no issues here.

Regarding the other discussion thread, it seems strange that after the swap only the LJ energies go to infinity, and not the Coulomb ones; that might be a bug somewhere in the exchange code. Two more sanity checks I would try are (1) try running REST with a much larger exchange interval, e.g. 5000 steps, and (2) see if switching to a different combination of Gromacs/Plumed helps.

Also, for testing purposes it's best to limit yourself to 2 replicas to try to set up a minimal working system.

Hi Miłosz Wieczór ,

I reconstructed a testing system that contains only 2 replicas each at 300K and 310K. Again the system I'm trying to simulate is just lysozyme in H2O with OPLS-AA force field and the .gro file I used was energy minimized after a 360-ns long regular MD production run (should be quite stable). I have attached the .mdp file here where I reduced the output frequency to 1. The submission file is also attached and I'm using an exchange frequency of 5000 this time as suggested. Again, I have got an error occurred at step 5001 saying that the total potential is infinite (data.out). I have checked both .log files and saw that their LJ and electrostatic energies were fine before swapping (topol*.log). Since I have tried many possibilities, I am suspecting it's the error in the original code. I am wondering if you could provide the .tpr files for your system so that I can test them (I tried to generate them using your script but I don't have "amber99bsc1.ff" force field in my share library....). Thank you again so much for your help!

Also the Gromacs version I use is 2018.6 with Plumed 2.5.2. Me and my colleague have tried different combination but this one seems to be the one that works without obvious patching issue.

I wasn't sure if I can share my system (it's a collaborative project), but I set up a quick DNA-based one where one of the DNA strands is treated as the solute, and verified that it works with my setup (gmx 2016.4, plumed 2.4.0) - there are 3 exchanges for the first 22 attempts. Please check if you keep getting the error on your machine.

Hi Miłosz Wieczór ,

Just a quick question, when I read the .top files generated with your script, it seems to me that only the partial charge and the LJ sigma (why it's not sigma?) are being scaled. I didn't see dihedral being scaled. Is there a mistake (or is it because your script does not support it as mentioned earlier)? Thank you!

I followed a slightly different philosophy here - my script duplicates all types found in the solute and adds scaled force field parameters in the parameter section. For instance in the [ atomtypes ] section, atomtype CE now has a corresponding atomtype YCE with altered epsilon.

Hi Miłosz Wieczór ,

I think I know which part you are talking about... regarding to my question, for me it seems like it's the LJ sigma that is being scaled rather than the LJ epsilon value.

Thanks for your patience!

If you look up e.g. the bottom of the [ atomtypes ] section in rest-0.968-topol.top, you will see rescaled epsilons and unchanged sigmas. Where is this screenshot taken from?

Hi Miłosz Wieczór ,

The screenshot was from the .top file I generated using your script with my system. I will check it more carefully again. Thanks!

Also I have used exactly the same plumed.dat, rest0.tpr, and rest1.tpr files from the folder you provided and modified the submission file as attached since we are using different version of Gromacs/Plumed and different machine. The only difference between my test and your test now should be the submission script/machine/Gromacs and Plumed version and I still got the following error:

Fatal error:

5308 particles communicated to PME rank 0 are more than 2/3 times the cut-off

out of the domain decomposition cell of their charge group in dimension x.

This usually means that your system is not well equilibrated.

The output file and .log files are attached. Now it looks like a version issue. Could I confirm with you again which version of Gromacs/Plumed you are using? Since in the submission file you provided it is Gromacs 5.1.4-Plumed 2.2.3 instead of gmx 2016.4, plumed 2.4.0. Thank you so much!

Thanks, I'll have to test it with OPLS topologies then - I only used it for CHARMM and Amber so far. That surprised me because rest-0.968-topol.top was also generated with my script and the formatting was completely different.

As you say, this might be a version issue indeed: if you collect enough evidence for that, it would be worthwhile to report a bug. Also, in hurry I accidentally sent you two submission scripts - for the tests I actually used the other one, submit-multi.sh, which runs mdrun_mpi from gmx2016.4-plu240.

Hi Miłosz Wieczór ,

We have installed Gromacs 2016.4 patched with Plumed 2.4.0 and run your .tpr files and again we still got the decomposition error as attached. We read the .log files and found the only difference should be the gcc we used. You are using gcc 7.3.0 while we are using gcc 9.1.0. I don't know if this is the thing that causes issue. I have attached the submission script I used.

Just want to make sure, you ran the .tpr files in the folder you provided and there was no such error right (or the infinite potential error)? Thank you so much!

Yes, I didn't run it for long, but long enough to have 3 exchanges and there was no DD error.

What you attached seems like a different error, though - previously the error occurred at the first exchange attempt, and now it reports a missing atom in a DD cell. Did you encounter exactly this one previously?

Hi Miłosz Wieczór ,

We have tried playing with different submission commands. Sometimes we got the domain decomposition error and sometimes we got the infinite potential error, we also got this "...could only obtain xxx atoms" error sometimes. Now I'm not sure where we did things wrong since we are using the same version of Gromacs/Plumed and also the exactly same inputs.....

Hi Miłosz Wieczór ,

I have tested again with your .tpr files and got the error like this again:

Fatal error:

5267 particles communicated to PME rank 4 are more than 2/3 times the cut-off

out of the domain decomposition cell of their charge group in dimension x.

This usually means that your system is not well equilibrated.

Another kind of error I got is attached here and it is related to the triclinic box...

I wonder if you have seen these kinds of error before. It says the system is not well equilibrated which is strange because you can run the same .tpr files in your machine.... Thank you so much again! We really appreciate your help.

Just to double-check, I ran the same system on a different HPC cluster with a completely different configuration (Gromacs 5.1.4, Plumed 2.3.0) and still no sign of abnormalities after the first exchanges were completed. At this point I'd start being suspicious about the compilation itself, or something peculiar in your cluster's architecture. I assume you have an experienced IT/HPC person around, often they'd knew much more about these things that I could ever imagine.

Still, I've seen the "not well equilibrated" error hundreds of times before - usually for good and tractable reasons - but I'm seeing the "could only obtain..." and the infinite LJ energy ones for the first time. Anyway, the "dlopen failed" log message looks a bit nasty, too - it's a dynamic library loader error, so maybe it's the first failure in a chain of errors.

Hi Miłosz Wieczór ,

Thank you so much for the help! I have reached out to the technical person of our cluster and we are trying to solve the issue now. In addition, I have also tried running your .tpr files on another cluster and it worked (Gromacs 2016.5/Plumed2.4.2)! Just when I tried to use my own system it failed again with "one or more water molecules can not be settled" at the exchange step as attached. Since my system is just lysozyme in water and I have done adequate energy minimization and NVT/NPT equilibrium for 100ps and 200ps, respectively, I'm not sure why I still have this unequilibrated issue.

I have consulted our technical person and he said that the "dlopen failed" error is harmless. Thank you!

Glad to hear you got a minimal example to work! With time most MD codes (including Gromacs) grew so complex that these seemingly unimportant compilation flags, architecture details or library issues can end up breaking more complex systems once in a while, so it's good to test more than one compilation if you suspect code as the culprit. Again, the HPC people usually know much more about these issues than most MD people (including me).

AFAIK the "water molecules can not be settled" is a SETTLE algorithm error. Functionally, SETTLE is the equivalent of LINCS but used exclusively for water molecules, and here both errors seem to be a byproduct of the exchange-induced crash at t=500 steps.

Hi Miłosz Wieczór ,

Yes, you are absolutely right!

One more thing, could you let us know how you prepared the .gro file for the generation of .tpr files? like, did you do energy minimization, then NVT equilibration, then NPT equilibration, and then energy minimization? I just want to make sure that my system will be well equilibrated enough. Thank you!

The file was prepared with simple energy minimization. I don't normally do the NVT step, just go for NPT and remove the first several ns from analyses to let it relax. I think the only case where a little more care is required is membrane systems.