[gmx-users] number of degrees of frredom

[Justin A. Lemkul]

Moeed wrote:
> Dear Justin,
> 
> The actual box I got by replicating monomer units has size of 15 0.5 
> 0.5!!!! That is chain has no space to move. I know you have already 
> given me some tips on high repulsion in the system but that had to do 
> with incorrect topology file. I assure you I have learned a lot from 
> your tips ;)

I'm glad to hear it.  The same principal still applies, though.  If a large 
system is failing, break it down into manageable pieces to determine where 
things start to break down.  It is easier to solve one problem at a time than 
several.

> This time I need to have such packed system because that gives me the 
> density I want. 0.9g/cm3 for PE! for a chain of 60 ethylene units the 
> size should be identical to the actual size 15 0.5 0.5. And I am 
> wondering how can I work with much longer chains ( 20 times longer)...
> 

You've been given advice about this.  You probably will have to simulate within 
a larger box and slowly compress your system.  You have constructed a linear 
chain of PE within a pre-defined box that gives you the density you want.  Ask 
yourself this: is that what PE looks like, completely rigid and extended? 
Certainly not.  The better approach would be as Vitaly and Chris have suggested 
- take a single molecule of PE in a very large box (forget your density for a 
moment) and simulate it.  It will probably adopt a much more reasonable 
structure that you can then work with.  Then slowly build up your system - 2, 4, 
8 units - and it will likely be far more stable.  If your box fits your system 
exactly, your atoms will be extremely close to one another and you will get a 
huge vdW repulsion.  If you've learned anything from me, you should already know 
that I want you to break down your energy terms to see what's failing...

> md grompp is giving error about degrees of freedom. Any advise on this 
> in appreciated...
> 
> Thank you,
> 
> md 
> grompp****************************************************************************
> aking dummy/rest group for Acceleration containing 2896 elements
> Making dummy/rest group for Freeze containing 2896 elements
> Making dummy/rest group for Energy Mon. containing 2896 elements
> Making dummy/rest group for VCM containing 2896 elements
> Number of degrees of freedom in T-Coupling group Eth is 5373.22
> Number of degrees of freedom in T-Coupling group EthB is 211.89
> Number of degrees of freedom in T-Coupling group EthE is 211.89
> Making dummy/rest group for User1 containing 2896 elements
> Making dummy/rest group for User2 containing 2896 elements
> Making dummy/rest group for XTC containing 2896 elements
> Making dummy/rest group for Or. Res. Fit containing 2896 elements
> Making dummy/rest group for QMMM containing 2896 elements
> T-Coupling       has 3 element(s): Eth EthB EthE
> Energy Mon.      has 1 element(s): rest
> Acceleration     has 1 element(s): rest
> Freeze           has 1 element(s): rest
> User1            has 1 element(s): rest
> User2            has 1 element(s): rest
> VCM              has 1 element(s): rest
> XTC              has 1 element(s): rest
> Or. Res. Fit     has 1 element(s): rest
> QMMM             has 1 element(s): rest
> 

There is no error message here.  What's the problem?

> mdpfile**************************************************************************************
> 
> title               =  PE-Hexane
> ;define                   =  -DPOSRES   ; tells gromacs to perform 
> position restrained dynamics/include posre.itp into topology used for 
> position restraint
> pbc                       =  xyz                ; use priodic BCs in all 
> directions
> 
> ;               Run control
> integrator          =  md               ; type of dynamics algorithm. 
> Here md uses a leap-frog algorithm for integrating Newtons's eq of motion
> dt                  =  0.002            ; in ps !
> nsteps              =  5000             ; length of simulation= nsteps*dt
> nstcomm             =  1                ; frequency for center of mass 
> motion removal
> 
> ;               Output control
> nstenergy           =  100              ; frequency to write energies to 
> energy file. i.e., energies and other statistical data are stored every 
> 10 steps
> nstxout             =  100              ; frequency to write 
> coordinates/velocity/force to output trajectory file. how often 
> snapshots are collected= nstxout*dt
> nstvout             =  100
> nstfout             =  0
> nstlog              =  100              ; frequency to write energies to 
> log file
> nstxtcout                 =  10         ; frequency to write coordinates 
> to xtc trajectory
> 
> ;               Neighbor searching
> nstlist             =  10               ; frequency to update neighbor 
> list. Neighborlist will be updated at least every 10 steps. Manual p80
> ns_type             =  grid             ; make a grid in the box and 
> only check atoms in neighboring grid cells when constructing a new 
> neighbor list every nstlist steps
> 
> ;               Electrostatics/VdW
> coulombtype         =  Shift        ; tells gromacs how to model 
> electrostatics. Coulomb/LJ potential is decreased over the whole range 
> and forces decay smoothly to zero betw$
> vdw-type                =  Shift        ; rcoulomb-switch/rvw-switch & 
> rcoulomb/rvdw
> rcoulomb-switch     =  0                ; where to start switching the 
> Coulomb potential
> rvw-switch          =  0                ; where to start switching the 
> LJ potential
> 
> ;               Cut-offs
> rlist               =  1.1              ; in nm. Cut-off distance for 
> short-range neighbor list
> rcoulomb            =  1.0              ; distance for coulomb cut-off
> rvdw                =  1.0              ; distance for coulomb cut-off
> 
> ;               Temperature coupling
> Tcoupl              =  berendsen
> tc-grps             =  Eth  EthB  EthE  ;HEX            ; groups to 

OK, this is a terrible idea.  You've broken down a bonded molecule into pieces, 
something that is completely unphysical.  Imagine if we coupled individual amino 
acids of a protein separately...

I'll assume this is the source of your error message, but if that's not the 
case, please be sure to actually provide all the relevant information.

-Justin

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

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

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