OPES_EXPLORE + MTMB error

[Cesare Malosso]

Hi all I have a problem running a MTMB simulation with OPES. Here is the plumed.dat:

 # vim:ft=plumed

 # UNITS LENGTH=A TIME=ps ENERGY=eV

 energy: ENERGY

 vol: VOLUME

 ecv1: ECV_MULTITHERMAL_MULTIBARIC ...

   ARG=energy,vol

   TEMP=900

   MIN_TEMP=300

   MAX_TEMP=1400

   PRESSURE=0.06022140857*6e5     # 50GPa 

   MIN_PRESSURE=0.06022140857*3e5 # 30GPa  

   MAX_PRESSURE=0.06022140857*1e6 # 100GPa 

 ... 

 opes: OPES_EXPANDED ARG=ecv1.* PACE=200 WALKERS_MPI FILE=DeltaFs.data OBSERVATION_STEPS=250

 PRINT STRIDE=500 ARG=* FILE=COLVAR

while here is the lammps thermo output right after the 250 observation steps:

Step Time KinEng PotEng TotEng Temp Press Density Volume

2050000           10    31.056984   -21070.689   -21039.632    942.22464    555692.63    1.9055235    948.12074 

 2050050        10.01    26.273216   -21070.977   -21044.704    797.09195    747404.05    1.8762234    962.92708 

 2050100        10.02    25.737338   -21086.413   -21060.675    780.83417    767262.22     1.798008    1004.8155 

 2050150        10.03    77.633003 -5.5752844e+14 -5.5752844e+14    2355.2747    716295.46    1.6529701    1092.9819 

 2050200        10.04    59.313079 -6.2772122e+29 -6.2772122e+29    1799.4743   -228117.97    1.6633667    1086.1504 

 2050250        10.05    60.289629 -7.0675126e+44 -7.0675126e+44    1829.1015    326762.84    1.8941914    953.79291 

 2050300        10.06    53.565966 -7.9573118e+59 -7.9573118e+59    1625.1151    686617.54    2.0448898    883.50302 

 2050350        10.07    50.899873 -8.9591366e+74 -8.9591366e+74    1544.2296    715774.25     2.095232     862.2751 

 2050400        10.08     42.97404 -1.0087091e+90 -1.0087091e+90    1303.7711    677631.31    2.0888116    864.92545 

 2050450        10.09    37.397087 -1.6999513e+104 -1.6999513e+104    1134.5743    642343.33    2.0548278    879.23005 

 2050500         10.1    34.879343 -2.7559656e+118 -2.7559656e+118    1058.1896    588594.15    2.0179982    895.27648 

 2050550        10.11    30.542315 -4.4679786e+132 -4.4679786e+132    926.61032    490326.92    2.0003377    903.18069 

 2050600        10.12    29.624016 -7.2434987e+146 -7.2434987e+146    898.75043     554125.9    2.0120358    897.92954 

 2050650        10.13    29.092038 -1.1743179e+161 -1.1743179e+161      882.611    612175.21    2.0346091    887.96728 

 2050700        10.14    28.333003 -1.9038073e+175 -1.9038073e+175    859.58293    619194.38    2.0496601    881.44681 

 2050750        10.15    26.427965 -3.0864574e+189 -3.0864574e+189    801.78679    574618.95    2.0560864    878.69185 

 2050800        10.16     25.82736 -5.0037728e+203 -5.0037728e+203    783.56531    634235.83    2.0584159    877.69745 

 2050850        10.17    27.612181 -1.1247128e+217 -1.1247128e+217    837.71423    604670.39    2.0540465    879.56446 

 2050900        10.18    28.388733 -2.4281391e+230 -2.4281391e+230    861.27372    619614.26    2.0465889    882.76953 

 2050950        10.19    29.054894 -5.2421024e+243 -5.2421024e+243    881.48408    581859.73    2.0386924     886.1888 

 2051000         10.2    30.805027 -1.131716e+257 -1.131716e+257    934.58064    586001.33    2.0354766    887.58883 

 2051050        10.21    32.326711 -2.4432583e+270 -2.4432583e+270    980.74634    599968.78    2.0363388    887.21304 

 2051100        10.22    32.696162 -5.2747432e+283 -5.2747432e+283    991.95497    630897.14    2.0372196    886.82947 

 2051150        10.23    31.714551 -1.1387628e+297 -1.1387628e+297     962.1743    613438.87     2.035648    887.51411 

ERROR: Non-numeric pressure - simulation unstable (../fix_nh.cpp:1069)

Last command: run 10000000

where you can notice that the potential energy stars growing (negatively). Am I doing something wrong?

Thank you
Cesare

[Michele Invernizzi]

Hi Cesare,

What version of PLUMED are you using? The keywords seem to be from an old development version, in the stable PLUMED version they look like this https://www.plumed.org/doc-v2.8/user-doc/html/_e_c_v__m_u_l_t_i_t_h_e_r_m_a_l__m_u_l_t_i_b_a_r_i_c.html

The potential energy is expected to drop, but it should not go below the typical values at T=TEMP_MIN, see for example fig 4a of the PRX paper (it's not a multithermal simulation, but the behavior is similar) where it starts close to 12 then when the bias is switched on it jumps around zero and then quickly learns the correct bias and the values return up.

Assuming everything else is setup correctly (e.g. all walkers are well equilibrated in the same thermodynamic condition, the system is stable in the whole thermodynamic range, the thermostat and the barostat have conservative settings, etc), one thing you can try is to run the simulation at a temperature closer to your lower temperature, so that the initial jump will be smaller and won't knock out the simulation. For instance you can have a look at the temperature ladder generated in the DeltaFs.data file and pick the value in the middle.

Hope this helps.

Best,

Michele

--
You received this message because you are subscribed to the Google Groups "PLUMED users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to plumed-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/eb3e0cc2-0258-4876-92e9-da375ee17593n%40googlegroups.com.

[Cesare Malosso]

Hi Michele,

Thank you very much for your answer. I'm currently using PLUMED Version: 2.8.0-dev (git: f04c6219d). I think the problem was, as you said, the low temperature bound. I tried to split my MTMB simulation in two, one in the range [300-600K][30-100GPa] and another in the range [500-1400K][30-100GPa] and the simulations are running without problems. Probably the neural network potential I'm using didnt like such a jump in the temperature.
I take the chance to ask you another question about the time needed for the bias to converge. I'm checking the rct in the DeltaFs.data output file, and I was wondering which was the tipical time needed to the bias in order to converge. I've attached some of the plots of the rct I've obtained and still after 10 ns it is not a constant. 

Thank you in advance

Best,

Cesare

[Michele Invernizzi]

Hi Cesare,

I suggest if possible to switch to the stable v2.8 version (or the current master branch), since there are some changes that might be useful for you. In particular the new default is to use geometric spacing between the temperatures (see here), which can help having a better overlap and thus improve convergence. Using the latest stable version will also make your simulations more future-proof.

Regarding the instabilities, yes, neural network potentials are typically more fragile, especially when pushed to configurations never seen during training, e.g. by an abruptly changing bias. But I am happy to hear that you managed to make it work by splitting the simulation in two. 

Another possibility would be to give an initial guess to the DeltaFS, by creating a DELTAFS file (or a STATE file) with all the temperatures and the guessed values, and then start the simulation from there using the RESTART=YES keyword. This effectively puts a limit to the bias, so that at the start it can be more gentle. If you try this, remember that the guess should be higher than the real DeltaFs, otherwise you will likely never visit that region again.

About checking for convergence, my suggestion is also to keep an eye on the energy and the volume, they should visit the whole range, otherwise the bias is not yet a good one. Another thing I find very informative, is the effective sampling size (or sampling efficiency) as a function of temperature/pressure (see here for how to calculate it). This tells you if the sampling is going well. Sometimes even if the rct is not flat one can still get a good sampling efficiency for the reweighting simply by dropping an initial transient. IThe important thing for reweighting is that the bias changes slowly (adiabatically) not that it is perfectly static.

Let me know if you have further questions!

Best,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/37da7bbe-7057-45a4-a307-18a998151066n%40googlegroups.com.

[Cesare Malosso]

Hi Michele,

Thank you very much for your answer, I'll check the effective sampling size then. 

Actually I've just found out that one of the two simulations (the one at [30-100GPa][300-600K]) is going "bad", I was checking how the system was exploring different volumes and I saw something a bit strange (attached plot "Volume vs Time"), and visualizing the trajectory it seems like atoms started drifting (attached a small movie of the system that is solid/plastic NH3). Is this what it should be expected?

Thank you.

Best,

Cesare

[Michele Invernizzi]

Hi Cesare,

It looks like you are not fixing the center of mass in your simulation. If you are using LAMMPS this is the relevant fix https://github.com/invemichele/OPESexpanded/blob/master/sodium/input.lmp#L52

Best,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/992276a8-12c4-4770-9f30-6f1a6914e4f9n%40googlegroups.com.

[Cesare Malosso]

Hi Michele,

yes, you're right, I thought It was enough to fix momentum zeroes at the beginning of the simulation, after the thermalisation (like I normally do in my NVT+NVE simulations).

Thank you

Best,

Cesare

[Cesare Malosso]

Hi Michele,

Sorry for bothering you again but I've tried to compute the effective sample size in my MTMB simulation (ranging from 80-100GPa, 500-1400K). I modified the weights since I'm running a NPT + MTMB in accordance with the PLUMED doc. Here you can take a look at the notebook with the results. The effective sampling size as a function of the temperature (here I've selected the target pressure equal to the one of the NPT simulation, p = 90GPa) doesn't look correct because the sampling is very poor in the target temperature interval, and it seems to sample also higher temperatures. Is there anything wrong?

Best,

Cesare

[Michele Invernizzi]

Hi Cesare,

There is a KBT missing after the bias when calculating the ESS, if you change that it will give you more reasonable results. 
This is with transient = len(colvar['opes.bias'])//10 #start

it's not very uniform, but you have more than 1% efficiency also at low temperature, so you should be able to get a nice reweighting (the statistical uncertainty goes like 1/sqrt(ESS) ). 

You can adapt this old script (https://github.com/invemichele/OPESexpanded/blob/master/chignolin/Analyze_neff-2D.py) to compute the ESS over both temperature and pressure.

It is probably not worth the time and the extra simulations, but if you want to get a more uniform ESS, you have to try to place more carefully the temperature/pressure steps. You can have a look at the replica exchange literature for the criteria to follow, several have been proposed. The geometric distribution over temperatures that I used assumes a constant heat capacity, which is often a bad approximation.

Best,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/65dad0d5-f7f3-426f-82db-7521a78d8e14n%40googlegroups.com.

[Cesare Malosso]

Hi Michele,

Thank you very much, I fixed the error on the weights and it should be fine. I also tried to compute the ESS on the entire range p-T thanks to the script you provided me and it doesn't look too bad for my MTMB-simulation in the range [500-1400K][80-100GPa]. 

I've also tried to run some MT-simulations and some simple NPT simulations in order to benchmark the MTMB simulation but I found some problems. Here I show you the average volume as a function of the temperature obtained from the MTMB simulation reweighted at a constant pressure of p=100GPa. I've also run a MT simulation at p = 100GPa and two NPT simulations at T=600K and T=900K, p = 100GPa. The results look good in this case, but when I turn to the average energy, there's no match between the results obtained from the different approaches. I've also found a mismatch between the energy predicted directly by LAMMPS and the energy computed by PLUMED. I've read that PLUMED does not include long tail corrections but it should not be the case since I didn't put them in the LAMMPS input file.

Thank you in advance,

Best,

Cesare

[Michele Invernizzi]

Dear Cesare,

I am happy to hear it's working nicely. 

The fact that the energies from PLUMED are different from those from LAMMPS may actually cause the problem and should be investigated further. It would be nice if you could share a minimum configuration that can be run and which presents this inconsistency.

Are the NPT averages obtained from the PLUMED energies as well? what happens if you redo all the ensemble averages and reweighting but using the energies from LAMMPS? Do you still get this inconsistency?
In the case of the missing tail corrections the difference is small enough that the reweighting can be done with the correct energies and the correct averages can be obtained without problems.

Best,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/de3dc686-294e-4fd5-ab6c-190ba381b4c7n%40googlegroups.com.

[Cesare Malosso]

Hi Michele,

Here you can find a small MTMB run of my system. File thermo1.out.0 contains the LAMMPS output quantities at the same step as in the plumed COLVAR file so you can directly compare the energies (eV in LAMMPS vs kJ/mol in PLUMED). 

NPT averages were obtained from the LAMMPS energies, just as time averages of the energy time series of the system.

Best,

Cesare

[Cesare Malosso]

Dear Michele,

Actually I found out that I was comparing the LAMMPS total energy of the system with the PLUMED potential energy, so of course there were different. Now, I'm considering the correct quantity (potential energy) via NPT simulation and I obtain something a bit more reasonable but still a bit annoying since it looks like the MT simulation is compatible with the NPT results for T < ~ 750K, while the MTMB for T> 750. Actually I was expecting the MTMB simulation to be less accurate at low temperature since the ESS is poor in that area, while I do not expect the MT to be not accurate at higher energy since the ESS looks pretty flat in the temperature interval. Also it is still a bit weird that the volume is compatible.

Best,

Cesare

[Michele Invernizzi]

Dear Cesare,

I did not have time yet to look at the files you shared, so it's nice to see that the problem of the energy mismatch is fixed.

Regarding the inconsistency shown in the figure, it's quite strange. Even if the adaptive bias is not 100% correct, the reweighting procedure should give the correct answer. I see possible explanations:

1. There is a bug somewhere in the reweighting script, and you are not really comparing the same quantities.
2. The uncertainty associated with these estimates are quite large, thus they are actually compatible. (you can estimate the uncertainty using block average or by running multiple independent simulations)
3. There is a subtle phase transition and the orange and blue lines are actually sampling two slightly different free energy basins (to check for this you should have a good look at the atomic configurations and check if something looks different between the two, and/or rerun the orange simulation starting from a blue configuration and vice versa.)

Hope this helps,

Michele

To view this discussion on the web visit https://groups.google.com/d/msgid/plumed-users/99d63a82-02c4-4755-ac0e-f22a3165bba5n%40googlegroups.com.