Bug #1399
B-state not set for position restraints when grompp -r filename is the same as the grompp -rb filename (or when neither -r nor -rb are provided)
Description
When doing slow-growth FEP (and perhaps other types of FEP), the position restraints are not read from the grompp -rb file if the same filename is passed to both the grompp -r and -rb options.
One might desire the same values for position restraints in the A and B states to keep the position restraint on and invariant of what is being changed in the FEP. I was doing this because I was inserting a ligand to the cognate ligand binding site of a protein and wanted to keep an allosteric binding site from changing conformation.
I used .mdp options like this:
free-energy = yes
init-lambda = 0
couple-lambda1 = vdw
couple-lambda0 = none
couple-intramol = no
couple-moltype = XOL
delta-lambda=0.0002
sc-alpha=0.5
sc-power=1
sc-r-power=6
sc-sigma=0.3
However, the restrained atoms were all pulled to the Cartesian origin as the simulation approached the B-state. This is because the posres_comB array was all zeroes. That's due to the following code in kernel/grompp.c:
static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi,
char *fnA, char *fnB,
int rc_scaling, int ePBC,
rvec com, rvec comB,
warninp_t wi)
{
int i, j;
read_posres (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
if (strcmp(fnA, fnB) != 0)
{
read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi);
}
}
I'll also note that if one does not define either -r or -rb, then again the posres_com array is fine (set to the -c positions) but the posres_comB values are all zeroes.
I'm not sure what the original purpose was with the above code, so I won't suggest any particular changes, but I don't think that the current behaviour is optimal, especially since it leads to -rb being ignored if the -r file has the same name as the -rb file, which seems like a common enough desire.
Now that I know what is going on, I am circumventing this by copying a.gro to b.gro and using grompp -r a.gro -rb b.gro
Chris.
Associated revisions
Always set b-state posres, even if identical to a-state
Previously we did not read the b-state position restraint
file (option -rb) if the name was identical to the a-state,
which caused the position restraints to be zero. Rather
than trying to be smart, we now always read it.
Fixes #1399.
Change-Id: I13d93ab667734f022e48143e4cd0672b7f303e1c
History
#1 Updated by Chris Neale about 7 years ago
- File example.tgz example.tgz added
Here is a tarball with test files. Read and run "run.sh"
#2 Updated by Erik Lindahl over 6 years ago
- Priority changed from Normal to High
#3 Updated by Erik Lindahl over 6 years ago
- Target version set to 5.0
#4 Updated by Gerrit Code Review Bot over 6 years ago
Gerrit received a related patchset '1' for Issue #1399.
Uploader: Erik Lindahl (erik@kth.se)
Change-Id: I13d93ab667734f022e48143e4cd0672b7f303e1c
Gerrit URL: https://gerrit.gromacs.org/3641
#5 Updated by Erik Lindahl over 6 years ago
- Status changed from New to Fix uploaded
#6 Updated by Erik Lindahl over 6 years ago
- Status changed from Fix uploaded to Resolved
#7 Updated by Erik Lindahl over 6 years ago
- Status changed from Resolved to Closed
Fix position restraint B-state fields
We now fill the B-state fields for position
restraints in the tpr, so test cases that use
position restraints must have their .tpr updated.
Refs #1399
Change-Id: Ic84f9d9627a59d9e32ec3ba817a345d5c76845f7