[gmx-users] questions regarding cutoff, charge groups, thermocouple and equlibration

[Justin A. Lemkul]

Kavyashree M wrote:
> Dear gromacs users,
> 1.  When one is using OPLSAA force field for simulation a protein of say 
> 100aa to 250aa long protein, TIP4P water model in gromacs 4.5.3
>     what are the values of cut offs to be used ie., rcoloumb, rvdw, 
> rvdw-switch and rlist? I have gone through the manual and papers and
>     tried using the following values:
> 
>    vdw_type = switch       | vdw_type = switch       | vdw_type = 
> switch           
>    coulomb_type = PME  | coulomb_type = PME  | coulomb_type = PME      
>    rlist = 1.2                     | rlist = 1.3                     | 
> rlist = 1.4                         
>    rcoulomb = 1.2            | rcoulomb = 1.3            | rcoulomb = 
> 1.4                
>    rvdw = 1.0                   | rvdw = 1.0                   | rvdw = 
> 1.0                      
>    rvdw_switch = 0.9        | rvdw_switch = 0.9        | rvdw_switch = 
> 0.9          
> 
>    but only last set of values did not give rise to any messages during 
> grompp regarding charge group radii being larger etc..  Are these 
> parameters acceptable?
>    I have gone through OPLS papers also but was unable to get a precise 
> answer for my question.
> 

It is difficult to say for certain how OPLS should be used with MD, since it was 
originally designed for different software doing MC simulations.  I commonly see 
rlist=rcoulomb=rvdw=1.0 with vdwtype = cutoff and coulombtype = PME in the 
literature.

> 2.  What is the type of thermocouple suggested during energy 
> minimization?  Berendson or Nose-hoover?
> 

During EM there are no temperatures.  Please consult some textbook material if 
this concept is foreign to you.

> 3. What is the way to confirm that a system has equilibrated after 
> starting the production MD? Is it RMSD from the initial structure, or 
> average structure?
>    I am attaching two rmsd plots of backbone after running for 12ns, 
> Does any of them indicates equilibration? 1st with 1.2nm (rlist = 
> rcoulomb=1.2nm);
>    2nd with 1.4nm, 1.2nm except 1.4nm gave the message regarding the sum 
> of charge group radii being larger
>    than the rlist - rvdw.
> 

There is no single metric that, in isolation, will tell you if a simulation is 
equilibrated.  You have to look at many variables and judge for yourself. 
First, is the ensemble stable (temperature, pressure, energy, whatever else)? 
Then, is the protein structure stable (RMSD, radius of gyration, secondary 
structure, hydrogen bonds, many more)?  From your RMSD plots, I would not be 
convinced that either are fully equilibrated, as one is still trending upwards. 
  More generally, as I said, you can't simply judge the simulation based on one 
criterion.

> 4. If MD is run on a protein which had intrinsic domain motion, how will 
> that equilibrate during MD? i.e., how will the rmsd plot appear when it 
> is equilibrated?

Depends entirely upon the types of motions that are involved.  Domain motions 
are usually very slow, and you may not see such movement on the nanosecond time 
scale with conventional MD.  Normal modes calculations are generally 
better-suited for these types of studies.

>     Will it equlibrate at all? will not the domain movement or movement 
> of a large loop in the structure prevent it from equilibration?
> 

Ideally, in the limit of infinite sampling, you would see your protein's 
structure interconvert between whatever conformations are expected at some 
predictable time interval.  I doubt you'll ever be able to run a conventional 
simulation long enough to see such behavior for even the smallest protein 
domains, unless you spend several years collected many microseconds of data. 
You may be able to see such motion faster with implicit solvent, but then the 
kinetics are meaningless in an absolute sense.

-Justin

> Kindly Give some suggestions and comments
> 
> Thanking you
> With Regards
> M. Kavyashree
> ------------------------------------------------------------------------
> 
> 
> ------------------------------------------------------------------------
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================