[gmx-users] Fwd: FEP calculations problem and posible BUG on GROMACS-3.3.2
Daniel Adriano Silva M
dadriano at gmail.com
Mon Oct 15 22:16:26 CEST 2007
Hello,
I have three problems here. First, I'm trying to made FEP calculations
in order to calculate the DG of binding for a small ligand (aminoacid)
to a protein, but my ligand (arginine) have a net charge of +1e, so
when I "disappear" this ligand the system takes a total charge of -1e
and I get the corresponding "note":
* State B has non-zero total charge: -1.000001e+00
My question is: what can I do?, because if at same time I try to
appear a Na+ atom I will be calculating the DG of binding of my ligand
plus the solvatation energy of one Na+ atom. In the other case, if I
doesn't compensate the charge I think PME will not work correctly.
A second problem is that I'm getting this:
WARNING 1 [file "topol_S.itp", line 39]:
Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to B
so I inferred that the problem is that I need also to specify a
B-state for bonds so for example I taked the line:
*; ai aj funct c0 c1 c2 c3
* 1 2 2 gb_2
and I transformed this to:
*; ai aj funct c0 c1 c2 c3
* 1 2 2 gb_2 gb_2
I made this for bonds, angles and dihedrals. My question is: is
correct my method?, and what is the meaning of this B-states?, it's to
say in what case I will want different bonds, angles, etc.. on my
simulation, and what happens If I doesn't specify this, will GROMACS
made the duplication of A state automatically?
Also I think it's important to mention that I have the topology of the
ligand: topol_S.itp, and the corresponding to the protein:
topol_E.itp. But I only get the warning for topol_S (disapearing
ligand), but not for the protein where that parameters on B-state are
not explicitly annotated. Question: why?
The third problem is what I think can be a BUG on 3.3.2. This warning
doesn't appear on 3.3.1:
*WARNING 2 [file aminoacids.dat, line 1]:
* T-Coupling group protein has fewer than 10% of the atoms (2381 out of
* 36249)
* Maybe you want to try Protein and Non-Protein instead?
I strongly think this is an error because in my MDP file I have explicitly:
* tc-grps = protein non-protein
and because in the GROMPP out also appears:
* T-Coupling has 2 element(s): Protein Non-Protein
I paste the output of grompp and my mdp:
GROMPP
#################################################################
creating statusfile for 12 nodes...
calling /usr/bin/cpp...
processing topology...
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
turning all bonds into constraints...
# G96ANGLES: 3528
# PDIHS: 1269
# IDIHS: 1158
# LJ14: 3898
# CONSTR: 2413
# SETTLE: 11289
Analysing residue names:
There are: 11290 OTHER residues
There are: 239 PROTEIN residues
There are: 0 DNA residues
Analysing Protein...
Analysing Other...
e frame at or first after this time.
-np int 12 Generate statusfile for # nodes
-[no]shuffle bool no Shuffle molecules over nodes
-[no]sort bool no Sort molecules according to X coordinate
-[no]rmvsbds bool yes Remove constant bonded interactions with virtual
sites
-load string Releative load capacity of each node on a
parallel machine. Be sure to use quotes around
the string, which should contain a number for
each node
-maxwarn int 10 Number of warnings after which input processing
stops
-[no]check14 bool no Remove 1-4 interactions without Van der Waals
-[no]zero bool no Set parameters for bonded interactions without
defaults to zero instead of generating an error
-[no]renum bool yes Renumber atomtypes and minimize number of
atomtypes
Back Off! I just backed up LAMBDA_0.0/mdout_lambda_0.0.mdp to
LAMBDA_0.0/#mdout_lambda_0.0.mdp.19#
checking input for internal consistency...
Generated 165 of the 1596 non-bonded parameter combinations
Excluding 3 bonded neighbours for Protein_E 1
Excluding 3 bonded neighbours for Protein_S 1
Excluding 2 bonded neighbours for SOL 11289
Excluding 1 bonded neighbours for NA+ 1
Excluding 1 bonded neighbours for CL- 0
NOTE:
State B has non-zero total charge: -1.000001e+00
processing coordinates...
double-checking input for internal consistency...
Velocities were taken from a Maxwell distribution at 310 K
renumbering atomtypes...
converting bonded parameters...
Walking down the molecule graph to make shake-blocks
initialising group options...
processing index file...
Opening library file
/global/home/est/dasm/PROGRAMAS/bin/../gromacs-3.3.2/build/share/gromacs/top/aminoacids.dat
WARNING 1 [file aminoacids.dat, line 1]:
T-Coupling group protein has fewer than 10% of the atoms (2381 out of
36249)
Maybe you want to try Protein and Non-Protein instead?
Making dummy/rest group for Acceleration containing 36249 elements
Making dummy/rest group for Freeze containing 36249 elements
Making dummy/rest group for VCM containing 36249 elements
Number of degrees of freedom in T-Coupling group Protein is 4729.80
Number of degrees of freedom in T-Coupling group Non-Protein is 67734.20
Making dummy/rest group for User1 containing 36249 elements
Making dummy/rest group for User2 containing 36249 elements
Making dummy/rest group for XTC containing 36249 elements
Making dummy/rest group for Or. Res. Fit containing 36249 elements
Making dummy/rest group for QMMM containing 36249 elements
T-Coupling has 2 element(s): Protein Non-Protein
Energy Mon. has 2 element(s): Protein Non-Protein
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
Checking consistency between energy and charge groups...
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 72x72x70, spacing 0.114 0.114 0.117
splitting topology...
Walking down the molecule graph to make shake-blocks
There are 12315 charge group borders and 11292 shake borders
There are 11292 total borders
Division over nodes in atoms:
3020 3021 3021 3021 3021 3021 3021 3021
3021 3021 3021 3019
writing run input file...
Back Off! I just backed up LAMBDA_0.0/lambda_0.0.tpr to
LAMBDA_0.0/#lambda_0.0.tpr.19#
There was 1 warning
#################################################################
#################################################################
MDP
#################################################################
title = dsilva Pos-Res
cpp = /usr/bin/cpp
integrator = md
dt = 0.002
nsteps = 1500000
pbc = xyz
nstlist = 10
rlist = 1.0
ns_type = grid
coulombtype = pme
rcoulomb = 1.0
vdwtype = cut-off
rvdw = 1.0
tcoupl = Berendsen
tc-grps = protein non-protein
tau-t = 0.1 0.1
ref-t = 310 310
Pcoupl = Berendsen
pcoupltype = isotropic
tau-p = 1.0
ref-p = 1.0
compressibility = 4.5e-5
fourierspacing = 0.12
pme_order = 4
optimize_fft = yes
ewald_rtol = 1e-5
gen_vel = yes
gen_temp = 310
gen_seed = 173529
constraints = all-bonds
constraint_algorithm = lincs
lincs_order = 4
; Free energy control stuff
free_energy = yes
init_lambda = 0.0
delta-lambda = 0
sc_alpha = 0.5
sc_sigma = 0.3
sc_power = 1.0
DispCorr = EnerPres
nstxout = 1000
nstvout = 1000
nstfout = 0
nstlog = 40
nstenergy = 40
energygrps = Protein non-protein
#################################################################
Any observation and/or recommendation will be welcome.
Thanks in advance.
Daniel Silva
More information about the gromacs.org_gmx-users
mailing list