exchange issues

[haoyuanchen]

In a 2x8x8 TUU run of alanine dipeptide with 6.56 force constant, in all the exchange steps in U dimension, there were exactly 16 accepted exchanges. Even if I repeat the exchange calculation by executing the script, I get exact the same result. In principle there should be some randomness. It also seems that within each group, only the first replica got exchanged. In one of the groups, there were actually 2 exchanges accepted according to the printout but in the pairs_for_exchange file, only 1 of the 2 were there.

The following table is from the first exchange step in U dimension and is organized as: replica_id replica_d1_param replica_d2_param replica_d3_param [accepted exchange pairs]

0 45.0 300.0 45.0 16 90.0 300.0 45.0 32 135.0 300.0 45.0 48 180.0 300.0 45.0 64 225.0 300.0 45.0 80 270.0 300.0 45.0 96 315.0 300.0 45.0 112 360.0 300.0 45.0 [96, 0] 1 45.0 300.0 90.0 17 90.0 300.0 90.0 33 135.0 300.0 90.0 49 180.0 300.0 90.0 65 225.0 300.0 90.0 81 270.0 300.0 90.0 97 315.0 300.0 90.0 113 360.0 300.0 90.0 [113, 1] 2 45.0 300.0 135.0 18 90.0 300.0 135.0 34 135.0 300.0 135.0 50 180.0 300.0 135.0 66 225.0 300.0 135.0 82 270.0 300.0 135.0 98 315.0 300.0 135.0 114 360.0 300.0 135.0 [114, 2] 3 45.0 300.0 180.0 19 90.0 300.0 180.0 35 135.0 300.0 180.0 51 180.0 300.0 180.0 67 225.0 300.0 180.0 83 270.0 300.0 180.0 99 315.0 300.0 180.0 115 360.0 300.0 180.0 [115, 3] 4 45.0 300.0 225.0 20 90.0 300.0 225.0 36 135.0 300.0 225.0 52 180.0 300.0 225.0 68 225.0 300.0 225.0 84 270.0 300.0 225.0 100 315.0 300.0 225.0 116 360.0 300.0 225.0 [100, 4] 5 45.0 300.0 270.0 21 90.0 300.0 270.0 37 135.0 300.0 270.0 53 180.0 300.0 270.0 69 225.0 300.0 270.0 85 270.0 300.0 270.0 101 315.0 300.0 270.0 117 360.0 300.0 270.0 [101, 5, 37, 53] 6 45.0 300.0 315.0 22 90.0 300.0 315.0 38 135.0 300.0 315.0 54 180.0 300.0 315.0 70 225.0 300.0 315.0 86 270.0 300.0 315.0 102 315.0 300.0 315.0 118 360.0 300.0 315.0 [118, 6] 7 45.0 300.0 360.0 23 90.0 300.0 360.0 39 135.0 300.0 360.0 55 180.0 300.0 360.0 71 225.0 300.0 360.0 87 270.0 300.0 360.0 103 315.0 300.0 360.0 119 360.0 300.0 360.0 [103, 7] 8 45.0 301.0 45.0 24 90.0 301.0 45.0 40 135.0 301.0 45.0 56 180.0 301.0 45.0 72 225.0 301.0 45.0 88 270.0 301.0 45.0 104 315.0 301.0 45.0 120 360.0 301.0 45.0 [104, 8] 9 45.0 301.0 90.0 25 90.0 301.0 90.0 41 135.0 301.0 90.0 57 180.0 301.0 90.0 73 225.0 301.0 90.0 89 270.0 301.0 90.0 105 315.0 301.0 90.0 121 360.0 301.0 90.0 [121, 9] 10 45.0 301.0 135.0 26 90.0 301.0 135.0 42 135.0 301.0 135.0 58 180.0 301.0 135.0 74 225.0 301.0 135.0 90 270.0 301.0 135.0 106 315.0 301.0 135.0 122 360.0 301.0 135.0 [122, 10] 11 45.0 301.0 180.0 27 90.0 301.0 180.0 43 135.0 301.0 180.0 59 180.0 301.0 180.0 75 225.0 301.0 180.0 91 270.0 301.0 180.0 107 315.0 301.0 180.0 123 360.0 301.0 180.0 [123, 11] 12 45.0 301.0 225.0 28 90.0 301.0 225.0 44 135.0 301.0 225.0 60 180.0 301.0 225.0 76 225.0 301.0 225.0 92 270.0 301.0 225.0 108 315.0 301.0 225.0 124 360.0 301.0 225.0 [108, 12] 13 45.0 301.0 270.0 29 90.0 301.0 270.0 45 135.0 301.0 270.0 61 180.0 301.0 270.0 77 225.0 301.0 270.0 93 270.0 301.0 270.0 109 315.0 301.0 270.0 125 360.0 301.0 270.0 [125, 13] 14 45.0 301.0 315.0 30 90.0 301.0 315.0 46 135.0 301.0 315.0 62 180.0 301.0 315.0 78 225.0 301.0 315.0 94 270.0 301.0 315.0 110 315.0 301.0 315.0 126 360.0 301.0 315.0 [126, 14] 15 45.0 301.0 360.0 31 90.0 301.0 360.0 47 135.0 301.0 360.0 63 180.0 301.0 360.0 79 225.0 301.0 360.0 95 270.0 301.0 360.0 111 315.0 301.0 360.0 127 360.0 301.0 360.0 [111, 15]

[haoyuanchen]

I'll look in to the Gibbs exchange codes in repex.

[antonst]

there were actually 2 exchanges accepted according to the printout

Which printout you are referring to?

[haoyuanchen]

Which printout you are referring to?

In the table I posted. In one group the accepted exchanges are [101, 5, 37, 53], but in the pairs_for_exchange file there was only [101 5] but no [37 53].

[antonst]

So this was a print out you have inserted correct?

[haoyuanchen]

So this was a print out you have inserted correct?

Yes.

[haoyuanchen]

I've run the simulation with different force constants (0.65, 6.5, 65), in all of them, the number of accepted exchange in each step in U dimension is 16, which correspond to a 25% acceptance rate. This is strange, since significantly different force constants should lead to significantly different acceptance rates.

[antonst]

I think this issue should be fixed with latest commit to devel. Please let me know it this is not the case.

[haoyuanchen]

The issue that only one of multiple accepted exchanges got recorded is fixed, and I'm seeing much more exchanges than before. However, it seems that even with a large force constant 65.6, the acceptance rate is still very high and is not significantly different from the simulation with a small force constant 0.656.

[antonst]

OK, I will try to spot check some things... but do you have any intuition why this is happening? What do you think may be causing this?

[haoyuanchen]

I'll also try to check the energies and print out something from the script. I'm not able to make a guess at this moment/

[haoyuanchen]

I think I found some glitches in the exchange code--not the exchange calculation itself but the way the function was called. I'm currently making changes and testing. Will make a commit if it turns out to be the case.

[haoyuanchen]

The problem is partially fixed now with my latest two commits--at least I can see the number of accepted exchanges decreases as the force constant increases. There're still some other problems that make the exchange rate not quite right and I'll continue to look into it.

[antonst]

There're still some other problems that make the exchange rate not quite right and I'll continue to look into it.

Thanks for your commits! Can you please elaborate a bit more? Where is the inconsistency?

[haoyuanchen]

Thanks for your commits! Can you please elaborate a bit more? Where is the inconsistency?

One of them is, one replica got "exchanged" multiple times in one cycle. I'm still trying to figure out how exactly Gibbs sampling should be executed in a single cycle within one group, but the current way seems incorrect anyway, which is, every replica will attempt to exchange with all others in the Gibbs way once. With this, N(N-1) exchanges were actually attempted in a group of N in one cycle. It should be N(N-1)/2 I believe. I'm fixing it now, but it's not as straightforward as just changing the inner loop over j to over j<i.

The other main problem now is in some cases some elements in the matrix don't look right, which lead to super-big exchange probabilities (actually Boltzmann factors). I'm still investigating this but a quick fix might work too.

[haoyuanchen]

I have a question for you by the way, if I see one replica appeared multiple times in a single pairs_for_exchange file, then how does the program handle this? For example, if these following two lines exist in the file:

1 2 1 3

Then, does 1 change with 2 first and then change with 3 (which is effectively 2 exchanged with 3), or just 1 change with 3, or just 1 change with 2?

[antonst]

1 will exchange with 2 first and than with 3, as a result 2 would have parameter from 1, 1 would have a parameter from 3 and 3 would have a parameter from 2.

[haoyuanchen]

OK, that's clear. In the previous implementation, if we only have 2 replicas A and B, in an exchange step, A will first try to exchange with B, and then B will try to exchange with A. This will cause an artificial increase in acceptance ratio because there should be only one exchange attempt between 2 replicas. My newest commit fixed this, although I'm not sure I fixed it in the right way...

[antonst]

btw, pairs which are in pairs_for_exchange file will definitely perform an exchange. What do you consider an attempt and at which stage this should be performed?

[haoyuanchen]

I know that pairs in that file will definitely exchange. Exchange "attempts" will happen to all pairs in a group, and only accepted exchanges go into that file. This was performed in global_ex_calculator.py, by something like

for g in groups
for replica in g
do_exchange

which will only return the accepted exchange attempts.

[antonst]

Is recording of attempts of exchanges something worth implementing?

[haoyuanchen]

Maybe not right now. By the way, it seems that the "double-counting" problem is not fixed properly although the exchange acceptance ratio looks better now. I'll continue to work on it.

[antonst]

I will fix "double-counting" problem, don't worry about that.

[haoyuanchen]

I will fix "double-counting" problem, don't worry about that.

I guess we have to be careful with that. If we just exclude all replicas that have already exchanged (accepted) from the pool (which is what I just did), that doesn't fully resolve the double-counting issue because if in the loop of replica 1, it attempted to exchange with replica 2 but failed, in the next loop of replica 2, it's still gonna attempt to exchange with replica 1, which is double-counting. In the gibbs_exchange routine, one replica actually tries (it tries, but maybe won't got the chance to go through all of them) to visit all other replicas attempting for exchange. So, it's hard to really define double-counting here.

[haoyuanchen]

On the opposite, if we ban all replicas that have "attempted" exchange instead of "accepted" exchange, then we won't be able to cover all possible permutations although we "fixed" the double-counting problem.

[antonst]

I thought you are still seeing multiple exchanges per replica in pairs_for_exchange files. This appears not to be the case now. So we can conclude that ''double counting'' problem is fixed now?

[haoyuanchen]

We could say that, but as I mentioned above, I think there's a more correct way to do it.

[antonst]

So if 1 attempted to exchange with 2 regardless of the outcome 2 should not attempt to exchange with 1?

[antonst]

I mean that if we are having 4 replicas and a 4x4 matrix, 1 replica attempts to exchange with everyone, replica 2 doesn't attempt to exchange with 1, replica 3 doesn't attempt to exchange with 1 and 2, etc. does this make sense?

[haoyuanchen]

So if 1 attempted to exchange with 2 regardless of the outcome 2 should not attempt to exchange with 1?

In principle (maybe) yes, but I'm not sure how to justify this, because in the gibbs_exchange routine, one replica actually tries (it tries, but maybe won't got the chance to go through all of them) to visit all other replicas attempting for exchange. If you look into the code (mainly line 58-62 in global_ex_calculator.py but also dependent lines nearby), you'll find that it's going through an array and will exit when a certain value is smaller than 0. I don't know if it exited at the 3rd element of the array, the 1st and 2nd element and their corresponding replicas are considered "attempted" or not. I guess I need to discuss with Taisung.

[haoyuanchen]

The "double-counting" problem actually is not a problem--it's ok to perform multiple rounds of exchanges in one step. The problem is calculating exchange acceptance ratio. Previously I just counted the number of pairs in the pair_for_exchange.dat file and divided by the total number of replicas. That is not correct. The right way to do it is, just compare the state of each replica before and after exchange. If it's different, then that replica got exchanged; otherwise it's not. Code will be like this:

accept_count = 0
states_before = []  # element i is the state of replica i before exchange
states_after = [] # element i is the state of replica i after exchange
for i in range(N):  # N replicas
    if states_before[i] != states_after[i]:
        accept_count += 1
acceptance_ratio = accept_count / N * 100

I'll implement this and correct some previous codes. Will give an update after I get some results.

[haoyuanchen]

With the new codes and the analysis script, the exchange acceptance ration of TUU seems to be quite reasonable. I'll do more tests and test TSU then. @antonst @taisung

[ibethune]

@haoyuanchen can you let us know if this issue was resolved?