[gmx-users] Domain decomposition error in alchemical free energy perturbation MD

2011-10-14 Thread Xiaoxiao He
Dear all,

I'm doing an "slow-growth" alchemical free energy perturbation calculation
of the formation of a disulfide bridge between two Cysteines with Gromacs.

I've had tried different ways to combine the topology of both state A and
state B, and finally settled with the most direct way -- to "mutate" the
atoms that have different partial charges in the two states, and transform
the two hydrogen atoms into "dummy" atoms. As suggested by Gromacs manual
and some messages from the internet, I put explicitly the OPLS parameters
for the bonds, pairs, angles and dihedrals changed from state A to state B,
and the grompp didn't give a warning. But when I was testing the production
simulation on two processors, there was a warning of fatal error in the log
file:

=
Initializing Domain Decomposition on 2 nodes
Dynamic load balancing: auto
Will sort the charge groups at every domain (re)decomposition
Initial maximum inter charge-group distances:
two-body bonded interactions: 3.774 nm, LJC Pairs NB, atoms 12 205
  multi-body bonded interactions: 0.640 nm, Ryckaert-Bell., atoms 1013 417
Minimum cell size due to bonded interactions: 4.151 nm
Using 0 separate PME nodes
Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
Optimizing the DD grid for 2 cells with a minimum initial size of 5.189 nm
The maximum allowed number of cells is: X 0 Y 0 Z 0

---
Program mdrun, VERSION 4.5.4
Source code file: domdec.c, line: 6436

Fatal error:
There is no domain decomposition for 2 nodes that is compatible with the
given box and a minimum cell size of 5.18876 nm
Change the number of nodes or mdrun option -rdd or -dds
Look in the log file for details on the domain decomposition
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

=

The atoms indexed 12 and 205 are a CD1 in Leu1 and another CB in Ala12
respectively. I checked the topology and found no entries containing both
these atoms. What's weird is that there's no way for these two atoms to have
bonded interactions, but it says in the log that they have a two-body bonded
interaction with a distance of 3.774 nm, which I cannot understand.

Does anyone have an explanation what this could mean? Any suggestion is
appreciated!

Thanks in advance!


Xiaoxiao He
Oct. 15, 2011



Attached are the input  parameters for the production md:
=
-Input Parameters:
   integrator   = sd
   nsteps   = 200
   init_step= 0
   ns_type  = Grid
   nstlist  = 10
   ndelta   = 2
   nstcomm  = 10
   comm_mode= Linear
   nstlog   = 100
   nstxout  = 10
   nstvout  = 10
   nstfout  = 10
   nstcalcenergy= 10
   nstenergy= 100
   nstxtcout= 1000
   init_t   = 0
   delta_t  = 0.0005
   xtcprec  = 1000
   nkx  = 36
   nky  = 36
   nkz  = 36
   pme_order= 4
   ewald_rtol   = 1e-05
   ewald_geometry   = 0
   epsilon_surface  = 0
   optimize_fft = TRUE
   ePBC = xyz
   bPeriodicMols= FALSE
   bContinuation= FALSE
   bShakeSOR= FALSE
   etc  = No
   nsttcouple   = -1
   epc  = Berendsen
   epctype  = Isotropic
   nstpcouple   = 10
   tau_p= 1
   ref_p (3x3):
  ref_p[0]={ 1.0e+00,  0.0e+00,  0.0e+00}
  ref_p[1]={ 0.0e+00,  1.0e+00,  0.0e+00}
  ref_p[2]={ 0.0e+00,  0.0e+00,  1.0e+00}
   compress (3x3):
  compress[0]={ 4.5e-05,  0.0e+00,  0.0e+00}
  compress[1]={ 0.0e+00,  4.5e-05,  0.0e+00}
  compress[2]={ 0.0e+00,  0.0e+00,  4.5e-05}
   refcoord_scaling = No
   posres_com (3):
  posres_com[0]= 0.0e+00
  posres_com[1]= 0.0e+00
  posres_com[2]= 0.0e+00
   posres_comB (3):
  posres_comB[0]= 0.0e+00
  posres_comB[1]= 0.0e+00
  posres_comB[2]= 0.0e+00
   andersen_seed= 815131
   rlist= 1.3
   rlistlong= 1.3
   rtpi = 0.05
   coulombtype  = PME
   rcoulomb_switch  = 0
   rcoulomb = 1.3
   vdwtype  = Switch
   rvdw_switch  = 0.8
   rvdw = 0.9
   epsilon_r= 1
   epsilon_rf   = 1
   tabext   = 1
   implicit_solvent = No
   gb_algorithm

[gmx-users] Re: Domain decomposition error in alchemical free energy perturbation MD

2011-10-15 Thread Xiaoxiao He
Hi everyone,

With the help of one senior colleague, the problem has been solved.

It turned out that if the parameter "couple-moltype" is turned on in the
mdp, there will be always a warning of distant pair interaction which should
not appear in the system. Once turned off, the warning disappears as well.
It was my misunderstanding of the functions of these parameters involved
(couple-moltype, couple-lambda0, couple-lambda1, etc) from the very
beginning. For a solvation/binding free energy calculation, where the
molecule is being decoupled from the system, it is more convenient to
provide these parameters because it saves the efforts of putting parameters
for state B in the topology.  But in the case of alchemical FEP calculation,
where only a few atoms of a residue are mutated, the decoupling of the whole
molecule becomes unnecessary and problematic. The correct way is
to implement  the topology for both state A and state B, which removes the
incorrect interactions and  is more reasonable.

Any insights are welcome. Thanks for your attention!

Best regards,
Xiaoxiao


On Sat, Oct 15, 2011 at 2:11 PM, Xiaoxiao He  wrote:

> Dear all,
>
> I'm doing an "slow-growth" alchemical free energy perturbation calculation
> of the formation of a disulfide bridge between two Cysteines with Gromacs.
>
> I've had tried different ways to combine the topology of both state A and
> state B, and finally settled with the most direct way -- to "mutate" the
> atoms that have different partial charges in the two states, and transform
> the two hydrogen atoms into "dummy" atoms. As suggested by Gromacs manual
> and some messages from the internet, I put explicitly the OPLS parameters
> for the bonds, pairs, angles and dihedrals changed from state A to state B,
> and the grompp didn't give a warning. But when I was testing the production
> simulation on two processors, there was a warning of fatal error in the log
> file:
>
>
> =
> Initializing Domain Decomposition on 2 nodes
> Dynamic load balancing: auto
> Will sort the charge groups at every domain (re)decomposition
> Initial maximum inter charge-group distances:
> two-body bonded interactions: 3.774 nm, LJC Pairs NB, atoms 12 205
>   multi-body bonded interactions: 0.640 nm, Ryckaert-Bell., atoms 1013 417
> Minimum cell size due to bonded interactions: 4.151 nm
> Using 0 separate PME nodes
> Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25
> Optimizing the DD grid for 2 cells with a minimum initial size of 5.189 nm
> The maximum allowed number of cells is: X 0 Y 0 Z 0
>
> ---
> Program mdrun, VERSION 4.5.4
> Source code file: domdec.c, line: 6436
>
> Fatal error:
> There is no domain decomposition for 2 nodes that is compatible with the
> given box and a minimum cell size of 5.18876 nm
> Change the number of nodes or mdrun option -rdd or -dds
> Look in the log file for details on the domain decomposition
> For more information and tips for troubleshooting, please check the GROMACS
> website at http://www.gromacs.org/Documentation/Errors
>
>
> =
>
> The atoms indexed 12 and 205 are a CD1 in Leu1 and another CB in Ala12
> respectively. I checked the topology and found no entries containing both
> these atoms. What's weird is that there's no way for these two atoms to have
> bonded interactions, but it says in the log that they have a two-body bonded
> interaction with a distance of 3.774 nm, which I cannot understand.
>
> Does anyone have an explanation what this could mean? Any suggestion is
> appreciated!
>
> Thanks in advance!
>
>
> Xiaoxiao He
> Oct. 15, 2011
>
>
>
> Attached are the input  parameters for the production md:
> =
> -Input Parameters:
>integrator   = sd
>nsteps   = 200
>init_step= 0
>ns_type  = Grid
>nstlist  = 10
>ndelta   = 2
>nstcomm  = 10
>comm_mode= Linear
>nstlog   = 100
>nstxout  = 10
>nstvout  = 10
>nstfout  = 10
>nstcalcenergy= 10
>nstenergy= 100
>nstxtcout= 1000
>init_t   = 0
>delta_t  = 0.0005
>xtcprec  = 1000
>nkx  = 36
>nky  = 36
>nkz  = 36
>pme_order