Dear Amber community,
I guess that there is a problem in the behavior of the "-ref" command line option in sander (AmberTools15-17).
>From Amber16 manual, I read the following:
"refc input (optional) reference coords for position restraints; also used for targeted MD"
What I understand from this is that the coordinates present in "refc" should be used to set an atom mask when "ntr > 1".
If that is correct, then I guess that this option does not work properly.
I have experienced the problem when setting a distance based mask (e.g., ":150.OG<:20").
As I want this mask to be consistent between different runs starting from different restart files,
I would like to base it on one single reference coordinate file.
I realized that the number of atoms that matches the mask differs from one run to another,
although the "refc' file remains the same.
For example:
With "restart.incprd" and "ref.inpcrd" being two slightly different sets of coordinates of the same system (no PBC), running:
sander -O -i md.in -p parm.prmtop -c restart.incprd -ref ref.inpcrd
or:
sander -O -i md.in -p parm.prmtop -c ref.incprd -ref ref.inpcrd
returns two different values for the number of atoms in the mask.
Actually, the second command returns the expected number of atoms, as given by the "ambmask" command when applied on "ref.inpcrd".
I tracked down the problem to mdread2.F90.
The reference coordinate file (refc) is read before assigning the "restraintmask" atom mask and the corresponding coordinates
are stored in "x(lcrdr)", while the input coordinates are stored "x(lcrd)" earlier in the code.
Then the mask is read using the "atommask" function, but with "x(lcrd)" as an input and not "x(lcrdr)".
Thus, the mask is systematically based on the input coordinates and not on the reference ones.
Correcting "lcrd" to "lcrdr" in the call of "atommask" fixes the problem and returns the expected mask.
More generally, I also realized that other mask based functions use the input coordinates to select the atoms
(i.e., all the calls to "atommask" in mdread2.F90 are based on "x(lcrd)").
The refc file is actually read only if ntr > 1 and not for other types of constraints/restraints.
I would for example expect "bellymask" to work on reference coordinates instead of input ones, but "refc" is not read when "ibelly = 1".
I fixed that locally to get the behavior that I expect, but I might be missing some more global vision of the code to be sure that it is correct.
I would greatly appreciate your opinion on this matter, and I apologies if I just misinterpreted the manual.
Best regards,
Antoine
________________________________
-- Dr. Antoine MARION
Postdoctoral Research Scientist
Department of Life Sciences
Theoretical Chemical Biology and Protein Modelling Group
Technische Universität München
Emil-Erlenmeyer-Forum 8
D-85354 Freising-Weihenstephan
Germany
Phone: +49 (0)8161 71-2271
www.linkedin.com/in/antoine-marion-b466844b<
http://www.linkedin.com/in/antoine-marion-b466844b>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri May 19 2017 - 09:30:03 PDT
