TMMC - biasing function not working as expected

[bartoszmazurwro]

Hi,

I think that the implementation of biasing function is not working as expected. I did some test after discovering that the macrostate histogram is not flat at all. So below is basic tmmc input as from your example:

{
  "SimulationType" : "MonteCarloTransitionMatrix",
  "NumberOfCycles" : 100000,
  "NumberOfInitializationCycles" : 50000,
  "NumberOfEquilibrationCycles" : 0,
  "PrintEvery" : 500,

  "ForceField" : "../ff/",

  "Systems" : [
    {
      "Type" : "Framework",
      "Name" : "../ff/MOF-303-E4D4",
      "NumberOfUnitCells" : [3, 2, 3],
      "ExternalTemperature" : 298.0,
      "ExternalPressure" : 1700,
      "ChargeMethod" : "Ewald",
      "UseChargesFrom" : "CIF_File",
      "CutOff" : 12.8,
      "MacroStateUseBias" : true,
      "MacroStateMinimumNumberOfMolecules" : 0,
      "MacroStateMaximumNumberOfMolecules" : 34
    }
  ],

  "Components" : [
    {
      "Name" : "../ff/tip4p",
      "FugacityCoefficient" : 1.0,
      "TranslationProbability" : 1.0,
      "ReinsertionProbability" : 1.0,
      "RotationProbability" : 1.0,
      "SwapProbability" : 2.0,
      "CreateNumberOfMolecules" : 0
    }
  ]
}

and the resulting CM stats are:

N	CM[-1]	CM[0]	CM[1]	bias	lnpi	Forward_lnpi	Reverse_lnpi	histogram
0	0	0.9798638317672098	1.0201361682327903	110.09154929561329	-110.09154929561329	0	0	2
1	8.33886158284123e-06	40.66821608283711	3.3317755783013037	95.28598687383239	-95.28598687383239	-0.6732110639064248	14.80556242178091	44
2	1.4581865729514916	21.39749553115384	3.144317895894669	94.98576819153085	-94.98576819153085	12.22487815644791	15.105781104082439	26
3	1.1629900684278927	32.60437177947679	1.2326381520953094	93.69391402056682	-93.69391402056682	12.993281547442235	16.397635275046465	35
4	3.041823285388791e-08	1.9999999695817672	1	78.63326943839479	-78.63326943839479	13.05144392515958	31.458279857218503	3
5	0.04912851003652767	21.72771999059909	3.22315149936438	73.49969014197472	-73.49969014197472	30.359667568550393	36.591859153638566	25
6	1.456991897680902	40.31125157258416	1.2317565297349404	72.16338020996022	-72.16338020996022	34.54334293619685	37.92816908565307	43
7	1.0225722552291159	48.50404984571708	2.4733778990538093	71.78721665515812	-71.78721665515812	34.37541019357265	38.30433264045518	52
8	1.9408785036173766	130.19582722187627	2.8632942745063277	70.5907415159586	-70.5907415159586	35.25867370870797	39.500807779654686	135
9	2.929844242237539	121.37452517987691	3.6956305778855536	70.66696248624635	-70.66696248624635	35.64750580741739	39.42458680936694	128
10	2.038186310190209	107.41168592509209	3.5501277647177307	70.19651408331289	-70.19651408331289	35.87970774236045	39.8950352123004	113
11	3.1293659596753223	84.25025171759081	2.620382322733864	70.29793905403235	-70.29793905403235	36.43463098648908	39.79361024158094	90
12	0.5155578985064321	31.31659214646694	1.1678499550266215	69.67541526446993	-69.67541526446993	36.25712080308529	40.416134031143365	33
13	0.6153193938302746	36.26596610749466	1.1187144986750588	68.89355844679392	-68.89355844679392	37.07479088266868	41.19799084881937	38
14	0.3484242956357137	51.8782901359305	1.773285568433792	67.37564600133544	-67.37564600133544	37.67258494616853	42.71590329427785	54
15	0.0029196449646635733	13.940204491656186	1.0568758633791502	62.247452510030264	-62.247452510030264	39.299753326942536	47.84409678558303	15
16	0.44896055047849714	81.76283383175567	2.7882056177658363	59.656713941105245	-59.656713941105245	45.191363842068185	50.434835354508046	85
17	1.0025401920651529	20.910474718275673	1.0869850896591735	59.94100971172605	-59.94100971172605	47.017582339188586	50.15053958388724	23
18	0.2722606096750329	38.10681075178639	1.6209286385385777	58.003221034872645	-58.003221034872645	47.0984532590382	52.088328260740646	40
19	0.002217090135358515	12.897544814320003	1.1002380955446387	52.45848424797649	-52.45848424797649	48.88244802530113	57.6330650476368	14
20	0.021604643079789786	52.249302176788845	1.7290931801313725	47.17818389415568	-47.17818389415568	55.08953432489896	62.91336540145761	54
21	1.2180068418684273	40.32941475754074	1.4525784005908284	47.05558651384753	-47.05558651384753	59.471978452641686	63.03596278176576	43
22	0	1.94836841768999	0.051631582310009805	50.44344644494345	-50.44344644494345	59.648102850669844	59.648102850669844	2
23	1.6200612365477468e-10	20.466522366574747	1.5334776332632467	28.465785924528223	-28.465785924528223	55.99133393665429	81.62576337108507	22
24	0.04464490340888382	38.918563331272296	1.0367917653188246	24.331395681171205	-24.331395681171205	78.9622590367757	85.76015361444209	40
25	1.0541973213442057	63.658315088912495	1.2874875897433091	23.847268933573933	-23.847268933573933	82.10740526451167	86.24428036203936	66
26	0.0028298165316442084	17.126153962495714	2.871016220972641	18.92095528833514	-18.92095528833514	82.30731833448792	91.17059400727815	20
27	2.015556590106407	91.1552512543473	2.829192155546293	16.998568699993985	-16.998568699993985	89.22952778478535	93.09298059561931	96
28	2.0243382377362553	94.05040202724405	3.925259735019711	16.622998342942825	-16.622998342942825	89.56862361765175	93.46855095267047	100
29	2.028379478806071	173.0896483000939	3.881972221099961	15.380587387468996	-15.380587387468996	90.23081329025463	94.7109619081443	179
30	2.7154700833395142	345.4943688709777	4.7901610456827095	14.344126888910171	-14.344126888910171	90.879919431203	95.74742240670312	353
31	2.1622785999414984	211.87356231922266	1.9641590808358782	14.03991507912615	-14.03991507912615	91.44751838185559	96.05163421648714	216
32	0.9591740408427766	400.6701700412596	8.370655917897569	12.682289383094332	-12.682289383094332	91.351420013844	97.40925991251896	410
33	5.761196680098511	26232.853043947518	528.3857593726829	8.558278330656794	-8.558278330656794	93.5178349995918	101.5332709649565	26767
34	19.229613561038665	4968892.700914816	99913.06948870377	0.00019677992102925537	-0.00019677992102925537	97.68906566322951	110.09135251569226	5068826

you can see that most of the time simulation was stuck in last macrostate. To test probable solutions I've tried to run (i) much longer simulation (2e6 production cycles), (ii) more recent rescaling of WL ("RescaleWangLandauEvery": 10), and (iii) additional equilibration ("NumberOfEquilibrationCycles" : 50000). In all cases simulation results were similar. As you can see from the long simulation statistics, the acceptance rate of swap move was extremely low:

Monte-Carlo moves statistics
===============================================================================

Component 0 [../ff/tip4p]
    Translation all:            13571086
    Translation total:           4523870    4526086    4521130
    Translation constructed:     4523870    4526086    4521130
    Translation accepted:        2264899    2267953    2260924
    Translation fraction:       0.500655   0.501085   0.500079
    Translation max-change:     0.370002   0.355597   0.412830

    Rotation all:            13566883
    Rotation total:           4521368    4523669    4521846
    Rotation constructed:     4521368    4523669    4521846
    Rotation accepted:        2259439    2266684    2260390
    Rotation fraction:       0.499725   0.501072   0.499882
    Rotation max-change:     0.628360   0.515731   0.555165

    Reinsertion(CBMC) all:            13570922
    Reinsertion(CBMC) total:          13570922
    Reinsertion(CBMC) constructed:    13506180
    Reinsertion(CBMC) accepted:            755
    Reinsertion(CBMC) fraction:       0.000056
    Reinsertion(CBMC) max-change:     0.000000

    Swap Insertion (CBMC) all:            13568232
    Swap Insertion (CBMC) total:          13568232
    Swap Insertion (CBMC) constructed:    13502811
    Swap Insertion (CBMC) accepted:           3513
    Swap Insertion (CBMC) fraction:       0.000259
    Swap Insertion (CBMC) max-change:     0.000000

    Swap Deletion (CBMC) all:            13571135
    Swap Deletion (CBMC) total:          13571135
    Swap Deletion (CBMC) constructed:    13571135
    Swap Deletion (CBMC) accepted:           3513
    Swap Deletion (CBMC) fraction:       0.000259
    Swap Deletion (CBMC) max-change:     0.000000

Based on the acceptance criteria for swap moves:

https://github.com/iRASPA/RASPA3/blob/f2c71c56ea792639ee14625dbc28a3a6426b9e47/src/raspakit/mc_moves/component/deletion.cpp#L121C1-L150C7

I noticed that the solution uses biasTransitionMatrix * Pacc. I haven't seen this approach before - in this and this paper, the biased acceptance criteria is calculated as:

$$ acc(old → new) = \min\left( 1, \frac{\pi(new) \exp[\eta(New)]}{ \pi(old) \exp[\eta(Old)]} \right) $$

where $\pi(n)$ is the probability of observing the system in microstate $n$, and $\eta(n)$ is a biasing potential for state $n$.

I’m curious to understand the reasoning behind the current implementation. Could you please share some insights on how biasTransitionMatrix * Pacc was derived? Any additional context would be really helpful.

Best
Bartosz

[Zhaoli2042]

Hi,

I suspect your TMMC bias wasn't updated enough times.

By default, the update (adjustBias function) happens every 1 million steps (updateTMEvery variable), and you had 100,000 MC cycles, with 34 molecules, this equals 3,400,000 MC steps. So your TMMC bias update happens only 3 times.

I would suggest running this for much longer.

Let me know if this works. My email is zhaoli2023@u.northwestern.edu

[bartoszmazurwro]

Hi Zhao,

thanks for reply. I also ran the simulation for 2,000,000 cycles which gave 67,848,258 steps and results are quite similar:

# performed: 67848258 steps
# collection matrix updated: 67 times
# minimum microstate: 0
# maximum microstate: 34
# column 1: N
# column 2: CM[-1]
# column 3: CM[ 0]0n# column 4: CM[+1]
# column 5: bias
# column 6: lnpi
# column 7: forward lnpi
# column 8: reverse lnpi
# column 9: histogram
N CM[-1] CM[0] CM[1] bias lnpi Forward_lnpi Reverse_lnpi histogram
0 0 4.893042497104846 2.106957502895155 76.9212675391497 -76.9212675391497 0 0 12
1 0.6219949057626123 7.290236933312044 1.0877681609253442 75.44988476936636 -75.44988476936636 -1.2006651839216445 1.4713827697833348 9
2 0.5034838364301751 8.496516163569826 1 74.57419254560442 -74.57419254560442 -0.6417137692102701 2.3470749935452804 10
3 0.48385799736062357 20.17712666501891 2.339015337620467 73.01531961347239 -73.01531961347239 0.044489900551234474 3.905947925677304 23
4 0.883344592486937 27.0104882299229 1.106167177590163 71.809748051954 -71.809748051954 1.6201837546696405 5.111519487195698 29
5 0 4 1 75.0761428351183 -75.0761428351183 1.845124704031392 1.845124704031392 5
6 0.2477511938546815 39.19584421341246 0.5564045927328599 71.60137100644263 -71.60137100644263 0.23568679159729178 5.3198965327070615 40
7 1.2226160647319153 118.61631287383227 5.161071061435796 71.24918916618502 -71.24918916618502 1.0447575139987104 5.672078372964679 125
8 1.5141876383168493 54.30557950229814 1.1802328593850056 70.80818659046926 -70.80818659046926 2.484908763661865 6.11308094868044 57
9 0 7.841510366778073 1.158489633221926 74.68552610083303 -74.68552610083303 2.2357414383166594 2.2357414383166594 9
10 0.3921063368587301 8.547359304536629 0.060534358604641465 73.60218677595041 -73.60218677595041 0.18563397738177612 3.3190807631992847 9
11 0.29563409554407494 85.75947556074917 3.9448903437067204 72.88551309424038 -72.88551309424038 -1.6826879784620772 4.035754444909309 90
12 1.347437466264623 41.464777773395994 1.1877847603393727 72.52691658851234 -72.52691658851234 0.9083659322184978 4.394350950637362 44
13 0.5564321547319123 21.316222952531582 1.127344892736506 72.4173119502142 -72.4173119502142 0.7822513430792966 4.503955588935497 23
14 0.03487784605102943 35.175077099222285 0.7900450547266847 68.49351857720121 -68.49351857720121 1.4883265885200458 8.427748961948474 36
15 0.052586178262547925 47.249180390485996 1.6982334312514507 65.47558055695167 -65.47558055695167 4.608564719644673 11.44568698219803 49
16 0.1450395837056476 34.83714030754432 1.0178201087500378 63.32354478118709 -63.32354478118709 8.08345523676309 13.597722757962615 36
17 0.001082654920313587 23.939299632348295 1.0596177127313926 56.84218570784973 -56.84218570784973 10.031867011558049 20.07908183129997 25
18 0 2 1 60.00315333800983 -60.00315333800983 16.918114201139872 16.918114201139872 3
19 0.022915800686499715 12.51882600743589 0.45825819187761024 54.76088764917338 -54.76088764917338 15.819501912471763 22.160379889976326 13
20 0.3031331092163384 28.220255521057005 0.47661136972666107 53.54528042703845 -53.54528042703845 18.81510801653046 23.375987112111254 29
21 5.300151587927945e-05 14.851827997641486 1.1481190008426352 45.03585135015195 -45.03585135015195 19.267637423432642 31.885416188997752 16
22 0.19714737427845996 15.540108960017738 2.2627436657038027 43.15613962555784 -43.15613962555784 29.25095241855368 33.76512791359186 18
23 1.0824367904532095 65.63772791270958 3.2798352968372124 41.060652840061614 -41.060652840061614 31.691334244037318 35.86061469908809 70
24 1.5208260785210788 81.53207238612904 2.947101535349864 40.086261239042706 -40.086261239042706 32.7999126637951 36.835006300107 86
25 0.1804111559776249 50.411642414919946 1.4079464291024295 37.79602582671741 -37.79602582671741 33.46148116081596 39.12524171243229 52
26 0.2078213018355901 16.78455242858774 1.007626269576669 36.94368888357181 -36.94368888357181 35.51613020335041 39.977578655577894 18
27 0.019165164904659283 8.9703010413967 1.0105337936986398 33.56921723616662 -33.56921723616662 37.09480423427197 43.35205030298308 10
28 0.2590348423842931 14.597815300121809 1.1431498574938996 31.73794220731648 -31.73794220731648 41.059943909843234 45.183325331833224 16
29 0.2573039992463463 40.20228487193517 1.5404111288184859 29.281576809446562 -29.281576809446562 42.54452409473064 47.63969072970314 42
30 1.3769808443026086 260.42320055047196 6.199818605225385 27.316103408125493 -27.316103408125493 44.33407045898127 49.60516413102421 268
31 5.374480337199476 269.2703065020251 6.355213160775371 27.12587757278979 -27.12587757278979 45.838697184914686 49.79538996635991 281
32 6.3698402345619805 100040.82219426648 2021.807965504125 21.258447221325305 -21.258447221325305 46.00631074331991 55.6628203178244 102069
33 4.138444324903763 4704558.591991897 96219.26957914741 11.2418634601484 -11.2418634601484 51.7390851742487 65.6794040790013 4800782
34 17.581676736049474 65046706.36488737 1296424.054921914 1.311423857562867e-05 -1.311423857562867e-05 61.76930373373652 76.92125442491113 66343149

Macrostates lower than 30 were visited less than 100 times. If I'm not mistaken, updateTMEvery cannot be set in simulation.json?

If you would like to reproduce the simulation I attach input files:
ff.zip
long.zip

[bartoszmazurwro]

Also, I know that for initial bias WL method is much faster but then TMMC converges faster (ref), maybe implementation of WL/TMMC would be beneficial?

[Zhaoli2042]

Thanks, I can try reproducing the results and see what happens. TMMC performs well/badly depending on the system, and water makes things a bit more difficult.

WL/TMMC is definitely something we can look into. I know other codes like FEASST has this already.

[salam1416]

@bartoszmazurwro @Zhaoli2042 : Were you able to figure out why this issue is happening and how to fix it?
I am facing similar issue with TMMC on pure water (TIP4P) and HKUST-1/MIL-47 with gRaspa.

[bartoszmazurwro]

I'm afraid this is due to poor biasing of pure TM. You could check out if this is also issue with WL/TM in FEASST software