About Brownian dynamics simulation

User discussions
1

GROMACS version:5.0.6
GROMACS modification: No
How does gromacs fix the temperature of the system in Brownian dynamics simulations? In the .mdp file, I have not provided any thermostats. From the source-code velocity should be in the order of 1/dt. So if I calculate temperature by using 0.5mv^2 = (3/2)kBT, then the temperature is not being maintained. Please help me to figure out this.

static void do_update_bd_Tconsts(double dt, real friction_coefficient,
int ngtc, const real ref_t[],
real rf)
{
/ This is separated from the update below, because it is single threaded */
int gt;

if (friction_coefficient != 0)
{
    for (gt = 0; gt < ngtc; gt++)
    {
        rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]/(friction_coefficient*dt));
    }
}
else
{
    for (gt = 0; gt < ngtc; gt++)
    {
        rf[gt] = sqrt(2.0*BOLTZ*ref_t[gt]);
    }
}

}

static void do_update_bd(int start, int nrend, double dt,
ivec nFreeze[],
real invmass[], unsigned short ptype[],
unsigned short cFREEZE[], unsigned short cTC[],
rvec x[], rvec xprime[], rvec v[],
rvec f[], real friction_coefficient,
real rf, gmx_int64_t step, int seed,
int gatindex)
{
/* note – these appear to be full step velocities . . . */
int gf = 0, gt = 0;
real vn;
real invfr = 0;
int n, d;

if (friction_coefficient != 0)
{
    invfr = 1.0/friction_coefficient;
}

for (n = start; (n < nrend); n++)
{
    real rnd[3];
    int  ng  = gatindex ? gatindex[n] : n;

    if (cFREEZE)
    {
        gf = cFREEZE[n];
    }
    if (cTC)
    {
        gt = cTC[n];
    }
    gmx_rng_cycle_3gaussian_table(step, ng, seed, RND_SEED_UPDATE, rnd);
    for (d = 0; (d < DIM); d++)
    {
        if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
        {
            if (friction_coefficient != 0)
            {
                vn = invfr*f[n][d] + rf[gt]*rnd[d];
            }
            else
            {
                /* NOTE: invmass = 2/(mass*friction_constant*dt) */
                vn = 0.5*invmass[n]*f[n][d]*dt
                    + sqrt(0.5*invmass[n])*rf[gt]*rnd[d];
            }

            v[n][d]      = vn;
            xprime[n][d] = x[n][d]+vn*dt;
        }
        else
        {
            v[n][d]      = 0.0;
            xprime[n][d] = x[n][d];
        }
    }
}

}

2

There is no velocity in Brownian dynamics.
The temperature is not maintained, but directly affect the magnitude of the random displacement.
GROMACS defines a “velocity” in BD as displacement/dt, but this is not really a velocity.

3

Thanks for the reply.
Can you kindly help me figure out how the temperature of the system reported in the log file is being calculated?

4

The 'kinetic energy" is computed as: 0.5 sum 0.5 gamma displacement^2 / dt.
This is not a real kinetic energy, as there are no velocities. But if the integrations is accurate, the computed temperature should be close to the requested temperature.

We should document this in the reference manual.

5

Thank you so much. It’s working perfectly.