Temperature fluctuation in brownian dynamics NVT simulation

[siladitya]

Hi,

I am running a NVT simulation with langevin thermostat. The problem I am facing is the temperature of the system during the course of the simulation is fluctuating a lot compare to the set temperature which should not be the case, as I am doing a NVT simulation. Earlier in my code I had a bead spring model without any harmonic angle potential with NVT ensemble and langevin thermostat and with that I never faced any problem regarding the temperature fluctuation. But recently I introduced a harmonic angle potential in my script with BD NVT ensemble and facing the problem of temperature fluctuation. I can't understand the problem and any help regarding this matter will be very helpful. I have also pasted my hoomd script.

from hoomd_script import *
import random
#from math import exp

#input data file
system = init.read_xml(filename="data.xml")

#define the bead-bead interactions

def LJ_tab(r, rmin, rmax, epsilon, sigma):
    V = 4 * epsilon * ( (sigma / r)**12 - (sigma / r)**6);
    F = 4 * epsilon / r * ( 12 * (sigma / r)**12 - 6 * (sigma / r)**6);
    return (V, F)

table = pair.table(width=1000)

table.pair_coeff.set('0', '0', func=LJ_tab, rmin=1.4288, rmax=5.358, coeff=dict(epsilon=46.9953, sigma=1.786))
table.pair_coeff.set('1', '1', func=LJ_tab, rmin=0.95, rmax=3.41, coeff=dict(epsilon=1.383656, sigma=1.1386458))
table.pair_coeff.set('2', '2', func=LJ_tab, rmin=0.95, rmax=3.24, coeff=dict(epsilon=2.947790, sigma=1.0809928))
table.pair_coeff.set('1', '2', func=LJ_tab, rmin=0.95, rmax=3.33, coeff=dict(epsilon=2.019586, sigma=1.1098193))
table.pair_coeff.set('0', '1', func=LJ_tab, rmin=1.16984, rmax=4.3869, coeff=dict(epsilon=8.06383, sigma=1.4623))
table.pair_coeff.set('0', '2', func=LJ_tab, rmin=1.14672, rmax=4.3002, coeff=dict(epsilon=11.77, sigma=1.4334))


#define the bond type and parameters
spring_constant_1 = 500
spring_constant_2 = 10000
ro_1 = 1.462
ro_2 = 0.444

bon = bond.harmonic()
bon.set_coeff('bond0', k=spring_constant_1, r0=ro_1)
bon.set_coeff('bond1', k=spring_constant_2, r0=ro_2)

ang = angle.harmonic()
ang.set_coeff('1-0-1', k=7000, t0=1.047197)

all = group.all()

# thermostat
#resize_time = 1e6
cool_time = 5e5
integrate.mode_standard(dt=0.001)
T_start = 75.0
T_end = 70.0
T_seed = 12345
bd1 = integrate.bdnvt(group=all, T=variant.linear_interp([(0, T_start), (cool_time, T_end)]), gamma_diam=False, seed=T_seed, limit=0.01)
bd1.set_gamma('0', gamma=1.786)
bd1.set_gamma('1', gamma=1.1386458)
bd1.set_gamma('2', gamma=1.0809928)

run(cool_time)

#box rescale
Lx2=31.201
Ly2=31.201
Lz2=31.201
resize_time = 5e5
box_resize = update.box_resize(Lx = variant.linear_interp([(0, 208.0), (resize_time, Lx2)]), Ly = variant.linear_interp([(0, 208.0), (resize_time, Ly2)]), Lz = variant.linear_interp([(0, 208.0), (resize_time, Lz2)]) )
box_resize.set_params(scale_particles = True)

# initialize particle velocities
random.seed(12345);
T = 70.0
for p in system.particles:
    mass = p.mass;
    vx = random.gauss(0, (T / mass)**(0.5))
    vy = random.gauss(0, (T / mass)**(0.5))
    vz = random.gauss(0, (T / mass)**(0.5))
    p.velocity = (vx, vy, vz)

update.zero_momentum(period=1e6)

# number of timesteps to run
run(resize_time)

#remove all the fixes so we don't accidentally redefine fixes
box_resize.disable()
bd1.disable()

#system equilibration
T_end = 70.0
T_seed = 12345
run_time = 5e5
#bd1 = integrate.bdnvt(group=all, T=variant.linear_interp([(0, T_start), (run_time, T_end)]), gamma_diam=False, seed=T_seed)
bd1 = integrate.bdnvt(group=all, T=T_end, gamma_diam=False, seed=T_seed)
bd1.set_gamma('0', gamma=1.786)
bd1.set_gamma('1', gamma=1.1386458)
bd1.set_gamma('2', gamma=1.0809928)

run(run_time)

#defining log file
analyze.log(filename='mylogsilicaT70_1.log', quantities=['potential_energy', 'temperature', 'kinetic_energy', 'pressure'], period=1000, overwrite=True)

# dump a .mol2 file for the structure information
mol2 = dump.mol2()
mol2.write(filename='tetheredsilica_particleT70_60deg_1.mol2')
dump.dcd(filename='tetheredsilica_particleT70_60deg_1.dcd', period=10000)

run(1e6)

#dump a xml file
xml = dump.xml()
xml.set_params(position=True, type=True, mass=True, velocity=True, bond=True, angle=True)
xml.write(filename="silicaparticle_temp70_60deg_1cooling.xml")


Thanks,
Siladitya

[Carolyn Phillips]

Here is a guess off the top of my head.

I wonder if your choice of harmonic spring angle coefficient is
introducing large velocity fluctuations, that the BD is constantly
"fighting" or "exciting".

I would try lowering the size of the time step of your simulation and
see what happens. In theory it is the tightest bond of your
simulation that determines your timestep size.

-Carolyn

> --
> You received this message because you are subscribed to the Google Groups
> "hoomd-users" group.
> To view this discussion on the web visit
> https://groups.google.com/d/msg/hoomd-users/-/yhaImFbd8AkJ.
> To post to this group, send email to hoomd...@googlegroups.com.
> To unsubscribe from this group, send email to
> hoomd-users...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/hoomd-users?hl=en.

[siladitya]

Thanks a lot for your suggestion, I will try that.

Regards,
Siladitya

[siladitya]

Hi Carolyn,

I tried lowering the size of the time step of my simualtion from dt=0.001 to dt=0.00001, but still the temperature of the system is not steady. Do you think I should try changing the spring constant (k) values.

Thanks,

Siladitya

On Wednesday, November 7, 2012 12:15:06 PM UTC-6, Carolyn wrote:

[Carolyn Phillips]

Yes, I would see what happens as a function of k and maybe look at the
coupling to the bath too.

> https://groups.google.com/d/msg/hoomd-users/-/4xejq45pNQ8J.

[Joshua Anderson]

One thing to be aware of is that NVT means constant AVERAGE T. The instantaneous kinetic energy T logged by hoomd *should* fluctuate. You can verify that it fluctuates correctly by testing scaling laws. If you double the system size, the std dev of T should decrease by 1/sqrt(2).
--------
Joshua A. Anderson, Ph.D.
Chemical Engineering Department, University of Michigan

[Eric Jankowski]

A bit off-topic, but at T=70, most of the LJ interactions you specify
will be negligible. Why not scale down the T's and epsilons so
they're closer to 1, and use something like Weeks-Chandler-Anderson to
model the excluded volume of the non-attractive beads? Also, if 7000
is really the sort of stiffness used to model your molecules (that is,
they are rigid about these triplets) you may want to consider using
rigid bodies instead. Hyper-stiff bonds and angles can lead to
pathological behavior; they are degrees of freedom that want to have
energy partitioned into them, but are too stiff to be able to do so
realistically.

Eric

[siladitya]

Thanks a lot for your suggestion.

Regards,
Siladitya