[gmx-users] au--protein separate equilibration

[Justin Lemkul]

On 1/17/13 12:15 PM, fatemeh ramezani wrote:
> 1. For separate equilibration of AU and Protein, I'm using this em.mdp file for equilibration in step 1, that Protein and SOL are frozen:

This is not an equilibration, this is an energy minimization, which will likely 
not accomplish much with so many particles frozen.

> title               =  n.pdb
> cpp                 =  /lib/cpp
> define              =  -DFLEXIBLE
> constraints         =  none
> integrator          =  steep
> dt                  =  0.002
> nsteps              =  40000
> ;constraint_algorithm  = shake
> ;shake_tol           = 0.0001
> ;nstenergy           =  10
> ;nstxtcout           =  10
> ;nstlist             =  5
> ns_type             =  grid
> rlist               =  1
> coulombtype         =  PME
> rcoulomb            =  1
> rvdw                =  1.2
> fourierspacing      =  0.12
> fourier_nx          =  0
> fourier_ny          =  0
> fourier_nz          =  0
> pme_order           =  6
> ewald_rtol          =  1e-5
> optimize_fft        =  yes
> ;       Energy minimizing stuff
> emtol               =  1000.0
> emstep              =  0.01
> freezegrps          =  Protein   SOL
> freezedim           =  Y Y Y    Y Y Y
> comm_mode           =  None
>

I suppose this .mdp file will allow whatever the Au particle is to move, but I 
don't see much point to running anything where the solvent is frozen.  Freezing 
water also makes -DFLEXBILE irrelevant.

> 2. And in step 2 of equilibration (that AU is freezed), em2.mdp file contain:
> title               =  n.pdb
> cpp                 =  /lib/cpp
> define              =  -DFLEXIBLE
> constraints         =  none
> integrator          =  steep
> dt                  =  0.002
> nsteps              =  40000
> ;constraint_algorithm  = shake
> ;shake_tol           = 0.0001
> ;nstenergy           =  10
> ;nstxtcout           =  10
> ;nstlist             =  5
> ns_type             =  grid
> rlist               =  1
> coulombtype         =  PME
> rcoulomb            =  1
> rvdw                =  1.2
> fourierspacing      =  0.12
> fourier_nx          =  0
> fourier_ny          =  0
> fourier_nz          =  0
> pme_order           =  6
> ewald_rtol          =  1e-5
> optimize_fft        =  yes
> ;       Energy minimizing stuff
> emtol               =  1000.0
> emstep              =  0.01
> freezegrps          =    AU   AUI     AUC
> freezedim           =  Y Y Y  Y Y Y  Y Y Y
> comm_mode           =  None
>
>
>
> 3. In simulation, I'm using this md.mdp file:
> title                    = Yo
> cpp                      = /lib/cpp
> define                   =
> integrator               = md
> tinit                    = 0.00000
> dt                       = 0.002
> nsteps                   = 250000
> nstxout                  = 2500
> nstvout                  = 2500
> nstlog                   = 500
> nstenergy                = 250
> nstxtcout                = 0
> pbc                      = xyz
> rlist                    = 1.00
> coulombtype              = PME
> r_coulomb                = 1.00

The correct keyword is "rcoulomb."

> fourierspacing           = 0.10

Why switch the fourierspacing here?  You used 0.12 before, and now you'r eusing 
0.10.

> pme_order                = 6
> ewald_rtol               = 1e-6
> vdw-type                 = switch
> rvdw-switch              = 0.90
> rvdw                     = 1.00

These vdw settings are all very different from the previous minimization step. 
It's usually a good idea to keep the settings the same throughout all steps, 
otherwise you're using some broken form of a force field intermittently.

> tcoupl                   = nose-hoover
> tc_grps                  = SOL AUS AUB AUI
> tau_t                    = 0.1 0.1 0.1 0.1
> ref_t                    = 300 300 300 300

Depending on the size of these groups, assigning them to separate thermostats 
may not be appropriate, and it doesn't appear that the protein is coupled 
anywhere.  Consult http://www.gromacs.org/Documentation/Terminology/Thermostats

> Pcoupl                   = no
> gen_vel                  = yes
> gen_temp                 = 10
> gen_seed                 = 173529
> ; GolP has been tested with lincs only

Then why turn off constraints?

> constraints              = none
> constraint_algorithm     = lincs
> lincs_order              = 4
> lincs_iter               = 1
> lincs_warnangle          = 80
> freezegrps               = AU AUI
> freezedim                = Y Y Y Y Y Y
> comm_mode                 = None
> Are these .mdp file right?
>

No.  See points above.  I can't comment on whether any of the settings are 
correct for use with the GolP force field, but there are numerous 
inconsistencies from a generic standpoint.

-Justin

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

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

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