[gmx-users] Re: free energy calculation , grompp crash

Justin A. Lemkul jalemkul at vt.edu
Mon Feb 7 19:23:19 CET 2011



Da-Wei Li wrote:
> hello
> 
> Here they are the command line and mdp file. I use Gromacs 4.5.3. This 
> is a test case only and the protein is 1UBQ. Grompp wills top for about 
> 10 minutes then go through.
> 

The efficiency of this kind of process will depend on the amount of available 
memory on the system.  You're asking grompp to decouple a huge amount of degrees 
of freedom, which will require a lot of memory to do.  It sounds like it's 
working, in any case, so there's no real problem.

Whether or not simultaneously decoupling the LJ and Coulombic interactions of a 
whole protein will generate a stable trajectory or sensible result is another 
matter.

-Justin

> ***********output of grompp*****************
> 
> Back Off! I just backed up mdout.mdp to ./#mdout.mdp.6#
> Generated 2278 of the 2278 non-bonded parameter combinations
> Generating 1-4 interactions: fudge = 0.5
> Generated 2278 of the 2278 1-4 parameter combinations
> Excluding 3 bonded neighbours molecule type 'Protein'
> turning H bonds into constraints...
> Excluding 2 bonded neighbours molecule type 'SOL'
> turning H bonds into constraints...
> Coupling 1 copies of molecule type 'Protein'
> Setting gen_seed to 8552
> Velocities were taken from a Maxwell distribution at 300 K
> ********************************
> 
> Command line and mdp file:
> 
> ******************************
> grompp -f pr1.mdp -c after_em.gro -t em.trr -p topol.top -o pr1
> ******************************
> define = -DPOSRES ; position restrain the protein
> ; Run parameters
> integrator = sd ; leap-frog integrator
> nsteps = 5000 ; 2 * 50000 = 100 ps
> dt = 0.002 ; 2 fs
> ; Output control
> nstxout = 1000 ; save coordinates every 2 ps
> nstvout = 5000 ; save velocities every 100ps 
> nstenergy = 1000 ; save energies every 2 ps
> nstlog = 1000 ; update log file every 2 ps
> ; Bond parameters
> continuation = no ; first dynamics run
> constraint_algorithm = lincs ; holonomic constraints 
> constraints = hbonds ; H bonds constrained
> lincs_iter = 1 ; accuracy of LINCS
> lincs_order = 4 ; also related to accuracy
> ; Neighborsearching
> ns_type = grid ; search neighboring grid cels
> nstlist = 10 ; 20 fs
> rlist = 0.8 ; short-range neighborlist cutoff (in nm)
> rcoulomb = 0.8 ; short-range electrostatic cutoff (in nm)
> rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
> ; Electrostatics
> coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
> pme_order = 4 ; cubic interpolation
> fourierspacing = 0.12 ; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl = V-rescale ; modified Berendsen thermostat
> tc-grps = Protein Non-Protein ; two coupling groups - more accurate
> tau_t = 0.1 0.1 ; time constant, in ps
> ref_t = 300 300 ; reference temperature, one for each group, in K
> ; Pressure coupling is off
> pcoupl = no ; no pressure coupling in NVT
> ; Periodic boundary conditions
> pbc = xyz ; 3-D PBC
> ; Dispersion correction
> DispCorr = EnerPres ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel = yes ; assign velocities from Maxwell distribution
> gen_temp = 300 ; temperature for Maxwell distribution
> gen_seed = -1 ; generate a random seed
> ;free energy stuff
> free_energy              = yes 
> init_lambda              = 0.0
> delta_lambda             = 0
> sc_alpha                 =0.5
> sc-power                 =1.0
> sc-sigma                 = 0.3 
> couple-moltype           = Protein 
> couple-lambda0           = vdw-q
> couple-lambda1           = none
> *******************************
> 
> thanks.
> 
> dawei
> 
> 
> 
> 
> On Mon, Feb 7, 2011 at 1:05 PM, TJ Mustard <mustardt at onid.orst.edu 
> <mailto:mustardt at onid.orst.edu>> wrote:
> 
>     Dawei,
> 
>      
> 
>     I have no problems with proteins in the thousands of atoms. Can you
>     post your command line and mdp files?
> 
>      
> 
>     Thank you,
> 
>     TJ Mustard 
> 
> 
>     On February 7, 2011 at 9:31 AM Da-Wei Li <lidawei at gmail.com
>     <mailto:lidawei at gmail.com>> wrote:
> 
>>     Well. It  actually isn't dead but becomes very slow for large
>>     proteins.   dawei
>>
>>     On Mon, Feb 7, 2011 at 11:44 AM, Da-Wei Li <lidawei at gmail.com
>>     <mailto:lidawei at gmail.com>> wrote:
>>
>>         hi,
>>         I did more test and found that it depended on size of the
>>         protein. Grompp will die when number of atoms of the protein
>>         is larger than about 200. Is it possible the source code limit
>>         the size of the protein that can be decoupled?
>>         thanks.
>>         dawei
>>
>>
>>         On Mon, Feb 7, 2011 at 10:34 AM, Da-Wei Li <lidawei at gmail.com
>>         <mailto:lidawei at gmail.com>> wrote:
>>
>>             Dear users
>>             I tried free energy calculation but grompp couldn't go
>>             through. It stops after
>>             *******************
>>             Generated 2278 of the 2278 non-bonded parameter combinations
>>             Generating 1-4 interactions: fudge = 0.5
>>             Generated 2278 of the 2278 1-4 parameter combinations
>>             Excluding 3 bonded neighbours molecule type 'Protein'
>>             turning H bonds into constraints...
>>             Excluding 2 bonded neighbours molecule type 'SOL'
>>             turning H bonds into constraints...
>>             Excluding 1 bonded neighbours molecule type 'CL'
>>             turning H bonds into constraints...
>>             Coupling 1 copies of molecule type 'Protein'
>>             *******************
>>             The CPU usage is 100%. 
>>             I just add following into the mdp file:
>>             ***************
>>             free_energy              = yes 
>>             init_lambda              = 0.0
>>             delta_lambda             = 0
>>             sc_alpha                 =0.5
>>             sc-power                 =1.0
>>             sc-sigma                 = 0.3 
>>             couple-moltype           = Protein  
>>             couple-lambda0           = vdw-q
>>             couple-lambda1           = none
>>             ***************
>>             Does anyone have some idea about this problem?  thanks.
>>             Another question is whether I can switch off "two
>>             molecules" (such as protein+ligand) in free energy
>>             calculation? I searched this list and got that 4.0.7 did
>>             support this. how about 4.5.4?
>>             dawei
>>
>>
>      
> 
>     TJ Mustard
>     Email: mustardt at onid.orst.edu <mailto:mustardt at onid.orst.edu>
> 
> 
>     --
>     gmx-users mailing list    gmx-users at gromacs.org
>     <mailto:gmx-users at gromacs.org>
>     http://lists.gromacs.org/mailman/listinfo/gmx-users
>     Please search the archive at
>     http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>     Please don't post (un)subscribe requests to the list. Use the
>     www interface or send it to gmx-users-request at gromacs.org
>     <mailto:gmx-users-request at gromacs.org>.
>     Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> 
> 

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

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

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



More information about the gromacs.org_gmx-users mailing list