[gmx-developers] Implementing Monte Carlo rejections in md.c

[j_haberstroh]

Hi Developers!

I'm working on including a MC step into my MD that conditionally accepts a
step with some probability that I compute from the difference between the
initial and final coordinates. This part is all fine and dandy, I'm getting
appropriate computations and appropriate probabilities.

I must do something beyond the simplest rejection of moves: I cannot simply
roll back x and v and expect my system to sample new positions; since the
eqns of motion are deterministic, I will never sample new positions, and the
MC will not produce any useful effect.

My approach right now is to save the positions and velocities into new
variables that I've snew'd called xumb and vumb, and then to write them back
to the state->x and state->v with copy_rvecn(). I do this MC step right
after the end of END UPDATE STEP 2. Then I want to scramble the velocities a
bit, probably by tcouple-ing.

My issue is that my system is totally indifferent to the value I set to
state->v. I've tried setting state->v to zero on rejection, tried
thermalizing at the beginning of the loop, at the end of the loop... nothing
I do seems to change the final position that I end up in.

Here's a little sample of some of my output:

dE_0=-57.107506...Step 29...dE_f=-57.140678...  Change in dE=-0.033173
dE_0=-57.107506...Step 30...dE_f=-57.140678...  Change in dE=-0.033173
dE_0=-57.107506...Step 31...dE_f=-57.140678...  Change in dE=-0.033173
dE_0=-57.107506...Step 32...dE_f=-57.140678...  Change in dE=-0.033173
dE_0=-57.107506...Step 33...dE_f=-57.140678...  Change in dE=-0.033173
Accepted!
dE_0=-57.140678...Step 34...dE_f=-57.187866...  Change in dE=-0.047188
Accepted!
dE_0=-57.187866...Step 35...dE_f=-57.242420...  Change in dE=-0.054554
Accepted!
dE_0=-57.242420...Step 36...dE_f=-57.296158...  Change in dE=-0.053738

Steps 29-33 are rejected over and over, and the final dE does not change
even though I have an additional update_tcouple() after each rejection.
Either something is conspiratorially abyssmal or the code is somehow
circumventing my attempts to brutally modify it.

Any tips? I want to pull a thermal ensemble of velocities.

Thanks!
John Haberstroh

--
View this message in context: http://gromacs.5086.x6.nabble.com/Implementing-Monte-Carlo-rejections-in-md-c-tp5014786.html
Sent from the GROMACS Developers Forum mailing list archive at Nabble.com.