[gmx-users] domain decomposition error in the energy minimization step
Kutzner, Carsten
ckutzne at gwdg.de
Fri Jan 13 09:55:31 CET 2017
Hi Qasim,
> On 12 Jan 2017, at 20:22, qasimpars at gmail.com wrote:
>
> Hi Carsten,
>
> I think I couldn't clearly explain the protocol that I follow. Sorry for that. Firstly, I do the EM, nvt (100 ps), npt (100 ps) and md (100 ns) steps for the equilibrium. In all those steps I use the below free energy parameters for the forward state:
>
> free-energy = yes
> init-lambda = 0
> delta-lambda = 0
> ...
>
> For the equilibrium simulations of the backward state I use the below free energy parameters in all MD steps:
>
> free-energy = yes
> init-lambda = 1
> delta-lambda = 0
> ...
>
> After that, I check whether the system has converged in both forward and backward states. If yes, 50 starting snapshots are extracted from the last 50 ns (converged part) of each of the two trajectories of 100 ns. Then 50 forward and 50 backward/reverse simulations of 50 ps each are carried out for the protein-ligand complex.
> For 50 ps forward simulations I use the above free energy parameters of the equilibrium forward simulation. As for 50 ps backward simulations I use its previous free energy parameters again.
>
> The answers to your questions:
> - In my snippets delta-lambda=0. Because I use 1-step switching process in all simulations. That is known as fast growth method also.
>
> If I change delta-lambda to 0.00002 or something like that, it would be slow growth method.
> -My system blows up in the md step of 100 ns backward state simulation with too many lincs error.
But I assume that in the 50 forward and backward simulations of 50 ps each that you spawn
from your equilibrated 50 ns trajectories you 'switch' lambda from lambda=0 at 0 ps to
lambda=1 at 50 ps, so you will have a delta-lambda != 0, right? Why else would you do the
switching simulations?
>
> Maybe the protocol I follow is wrong especally for the backward state?
Maybe the protein+ligand complex is not stable, at least not in your simulation?
For it to be stable, the ligand should stay in the binding pocket in an equilibrium
MD simulation, without any further constraints.
Carsten
> As far as I know, the slow growth method doesn't always provide the accurate results... I therefore prefer to use the fast growth method.
> Another interesting thing is that no one has any experience on this topic on the gmx user mailing list :(
>
> Any suggestions will be appreciated.
>
> Thanks in advance.
>
>> On 12 Jan 2017, at 18:34, "Kutzner, Carsten" <ckutzne at gwdg.de> wrote:
>>
>> Hi Qasim,
>>
>>> On 11 Jan 2017, at 20:29, Qasim Pars <qasimpars at gmail.com> wrote:
>>>
>>> Dear Carsten,
>>>
>>> Thanks. The forward state simulations works properly with mdrun -ntmpi 8
>>> -ntomp 2 or mdrun -ntmpi 4 -ntomp 4 as you suggested.
>>> For the backward state GROMACS still gives too many lincs warning error
>>> with those mdrun commands in the md step, indicating the system is far from
>>> equilibrium. I used the below free energy parameters with the em, nvt, npt
>>> and md steps for the backward state (fast growth). And to keep the ligand
>>> in the active site of protein I use some restraints under the
>>> intermolecular_interactions in the top file as I told in my previous email.
>>>
>>> free-energy = yes
>>> init-lambda = 1
>>> delta-lambda = 0
>>> sc-alpha = 0.3
>>> sc-power = 1
>>> sc-sigma = 0.25
>>> sc-coul = yes
>>> couple-moltype = ligand
>>> couple-intramol = no
>>> couple-lambda0 = vdw-q
>>> couple-lambda1 = none
>>> nstdhdl = 100
>>>
>>> My questions:
>>>
>>> 1) Do you think that I should use the above free energy paramaters in all
>>> MD steps (em, nvt, npt and md)?
>> I think you want them only in the MD part.
>>
>>>
>>> 2) The structure seems fine before arising the lincs warning. First linc
>>> warning belongs to protein atoms (not ligand atoms). However the ligand
>>> doesn't move out of the active site. Maybe the following mdp file for the
>>> md step is not correct or the backward state simulation of a complex
>>> structure is something impossible with GROMACS?
>> It should work in both directions.
>> Why is your delta-lambda zero in your snippets?
>>
>>>
>>> 3) Maybe the intermolecular_interactions section and couple- flags don't
>>> turn off intramolecular interactions of the ligand and turn on state B?
>>>
>>> 4) How can I get rid of the linc warnings in the md step?
>> How many LINCS warnings do you get?
>> Does you system blow up?
>>
>> Carsten
>>
>>>
>>> #md.mdp:
>>>
>>> ; RUN CONTROL
>>> integrator = sd
>>> nsteps = 50000000
>>> dt = 0.002
>>> comm-mode = Linear
>>> nstcomm = 100
>>>
>>> ; OUTPUT CONTROL
>>> nstxout-compressed = 500
>>> compressed-x-precision = 1000
>>> nstlog = 500
>>> nstenergy = 500
>>> nstcalcenergy = 100
>>>
>>> ; BONDS
>>> constraint_algorithm = lincs
>>> constraints = all-bonds
>>> lincs_iter = 1
>>> lincs_order = 4
>>> lincs-warnangle = 30
>>> continuation = no
>>>
>>> ; NEIGHBOUR SEARCHING
>>> cutoff-scheme = verlet
>>> ns-type = grid
>>> nstlist = 20
>>> pbc = xyz
>>>
>>> ; ELECTROSTATICS-EWALD
>>> coulombtype = PME
>>> rcoulomb = 1.0
>>> ewald_geometry = 3d
>>> pme-order = 4
>>> fourierspacing = 0.12
>>>
>>> ; VAN DER WAALS
>>> vdwtype = Cut-off
>>> rvdw = 1.0
>>> rvdw-switch = 0.9
>>> DispCorr = EnerPres
>>>
>>> ; TEMPERATURE COUPLING (SD - Langevin dynamics)
>>> tc_grps = Protein ligand
>>> tau_t = 1.0 1.0
>>> ref_t = 298.15 298.15
>>>
>>> ; PRESSURE COUPLING
>>> pcoupl = Parrinello-Rahman
>>> pcoupltype = isotropic
>>> tau_p = 2
>>> ref_p = 1.0
>>> compressibility = 4.5e-05
>>>
>>> ; VELOCITY GENERATION
>>> gen_vel = yes
>>> gen_seed = -1
>>> gen_temp = 298.15
>>>
>>> ;FREE ENERGY
>>> free-energy = yes
>>> init-lambda = 1
>>> delta-lambda = 0
>>> sc-alpha = 0.3
>>> sc-power = 1
>>> sc-sigma = 0.25
>>> sc-coul = yes
>>> couple-moltype = ligand
>>> couple-intramol = no
>>> couple-lambda0 = vdw-q
>>> couple-lambda1 = none
>>> nstdhdl = 100
>>>
>>> Thanks in advance.
>>>
>>>> On 11 January 2017 at 10:49, Kutzner, Carsten <ckutzne at gwdg.de> wrote:
>>>>
>>>> Dear Qasim,
>>>>
>>>> those kinds of domain decomposition 'errors' can happen when you
>>>> try to distibute an MD system among too many MPI ranks. There is
>>>> a minimum cell length for each domain decomposition cell in each
>>>> dimension, which depends on the chosen cutoff radii and possibly
>>>> other inter-atomic constraints. So this is normally just a technical
>>>> limitation and not a problem with the MD system.
>>>>
>>>> You can do the following steps to circumvent that issue:
>>>>
>>>> a) use less ranks (at the domain decomposition limit, the parallel
>>>> efficiency suffers anyway)
>>>>
>>>> b) use separate PME ranks, so that you get less and larger domains
>>>> on the MPI ranks that do domain decomposition
>>>> (use mdrun -nt 20 -npme 4 ... for example, you will have to try
>>>> around a bit with the exact number of PME ranks for optimal
>>>> performance - or use the gmx tune_pme tool for that!)
>>>>
>>>> c) In case you haven't done so already, compile GROMACS with
>>>> MPI *and* OpenMP. Then, by using MPI domain decomposition plus
>>>> OpenMP parallelism within each MPI rank, you will be
>>>> able to use more cores in parallel even for smaller MD systems.
>>>>
>>>> Use mdrun -ntmpi 10 -ntomp 2 for 10 ranks * 2 threads or
>>>> mdrun -ntmpi 4 -ntomp 5 for 4 ranks * 5 threads
>>>>
>>>> With real MPI, you would use something like
>>>>
>>>> mpirun -np 10 gmx mdrun -ntomp 2 ...
>>>>
>>>> Don't forget to check your simulation performance, there will
>>>> be better and worse choices in terms of these decomposition parameters.
>>>>
>>>> Happy simulating!
>>>> Carsten
>>>>
>>>>
>>>>> On 11 Jan 2017, at 08:33, Qasim Pars <qasimpars at gmail.com> wrote:
>>>>>
>>>>> Dear users,
>>>>>
>>>>> I am trying to simulate a protein-ligand system including ~20000 atoms
>>>> with
>>>>> waters using GROMACS-2016.1. The protocol I tried is forward state for
>>>> the
>>>>> free energy calculation. The best ligand pose used in the simulations was
>>>>> got by AutoDock. At the beginning of the simulation GROMACS suffers from
>>>>> domain decomposition error in the energy minimization step:
>>>>>
>>>>> Fatal error:
>>>>> There is no domain decomposition for 20 ranks that is compatible with the
>>>>> given box and a minimum cell size of 1.7353 nm
>>>>> Change the number of ranks or mdrun option -rdd
>>>>> Look in the log file for details on the domain decomposition
>>>>>
>>>>> I checked the complex structure visually. I don't see any distortion in
>>>> the
>>>>> structure. To check whether the problem is the number of nodes, I used 16
>>>>> nodes:
>>>>>
>>>>> gmx mdrun -v -deffnm em -nt 16
>>>>>
>>>>> The energy minimization step was completed successfully. For the NVT
>>>> step I
>>>>> used 16 nodes again.
>>>>>
>>>>> gmx mdrun -v -deffnm nvt -nt 16
>>>>>
>>>>> After ~2000000 steps the system gave too many lincs warning.
>>>>>
>>>>> Whereas there is no problem with wild type protein when it is simulated
>>>>> without using -nt 16. These domain decomposition error and lincs warning
>>>>> arise for complex structure.
>>>>>
>>>>> By the way to keep the ligand in the active site of protein I use the
>>>> bond,
>>>>> angle and dihedral restraints under [ intermolecular_interactions ]
>>>> section
>>>>> in the top file.
>>>>>
>>>>> [ intermolecular_interactions ]
>>>>> [ bonds ]
>>>>> 313 17 10 0.294 0.294 10.000 0.000 0.294
>>>>> 0.294 10.000 4184.000
>>>>>
>>>>> [ angle_restraints ]
>>>>> 312 313 17 313 1 140.445 0.000 1
>>>>> 140.445 41.840 1
>>>>> 313 17 19 17 1 107.175 0.000 1
>>>>> 107.175 41.840 1
>>>>>
>>>>> [ dihedral_restraints ]
>>>>> 300 312 313 17 1 56.245 0.000 0.000
>>>>> 56.245 0.000 41.840
>>>>> 312 313 17 19 1 -3.417 0.000 0.000
>>>>> -3.417 0.000 41.840
>>>>> 313 17 19 14 1 -110.822 0.000 0.000
>>>>> -110.822 0.000 41.840
>>>>>
>>>>> The mdp file used for the energy minimization is follows:
>>>>>
>>>>> define = -DFLEXIBLE
>>>>> integrator = steep
>>>>> nsteps = 50000
>>>>> emtol = 1000.0
>>>>> emstep = 0.01
>>>>> nstcomm = 100
>>>>>
>>>>> ; OUTPUT CONTROL
>>>>> nstxout-compressed = 500
>>>>> compressed-x-precision = 1000
>>>>> nstlog = 500
>>>>> nstenergy = 500
>>>>> nstcalcenergy = 100
>>>>>
>>>>> ; BONDS
>>>>> constraints = none
>>>>>
>>>>> ; NEIGHBOUR SEARCHING
>>>>>
>>>>> cutoff-scheme = verlet
>>>>> ns-type = grid
>>>>> nstlist = 10
>>>>> pbc = xyz
>>>>>
>>>>> ; ELECTROSTATICS & EWALD
>>>>> coulombtype = PME
>>>>> rcoulomb = 1.0
>>>>> ewald_geometry = 3d
>>>>> pme-order = 4
>>>>> fourierspacing = 0.12
>>>>>
>>>>> ; VAN DER WAALS
>>>>> vdwtype = Cut-off
>>>>> vdw_modifier = Potential-switch
>>>>> rvdw = 1.0
>>>>> rvdw-switch = 0.9
>>>>> DispCorr = EnerPres
>>>>>
>>>>> ; FREE ENERGY
>>>>> free-energy = yes
>>>>> init-lambda = 0
>>>>> delta-lambda = 0
>>>>> sc-alpha = 0.3
>>>>> sc-power = 1
>>>>> sc-sigma = 0.25
>>>>> sc-coul = yes
>>>>> couple-moltype = ligand
>>>>> couple-intramol = no
>>>>> couple-lambda0 = vdw-q
>>>>> couple-lambda1 = none
>>>>> nstdhdl = 100
>>>>>
>>>>> I removed the free energy lines in the em.mdp and [
>>>>> intermolecular_interactions ] section in the top file but GROMACS still
>>>>> gives the domain decomposition error for the complex structure.
>>>>>
>>>>> Will you please give suggestions on getting rid of the lincs warning and
>>>>> domain decomposition messages?
>>>>>
>>>>> I would appreciate any kind of help.
>>>>>
>>>>> Thanks.
>>>>>
>>>>> --
>>>>> Qasim Pars
>>>>> --
>>>>> Gromacs Users mailing list
>>>>>
>>>>> * Please search the archive at http://www.gromacs.org/
>>>> Support/Mailing_Lists/GMX-Users_List before posting!
>>>>>
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>
>>>>> * For (un)subscribe requests visit
>>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>>
>>>> --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at http://www.gromacs.org/
>>>> Support/Mailing_Lists/GMX-Users_List before posting!
>>>>
>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>
>>>> * For (un)subscribe requests visit
>>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>
>>>
>>>
>>> --
>>> Qasim Pars
>>> --
>>> Gromacs Users mailing list
>>>
>>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>>
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>
>>> * For (un)subscribe requests visit
>>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
>>
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list