Planarity restraint ignored when dist_esd equal to 0
Dear PHENIX developers, My colleague noticed that planarity restraint of a ligand was ignored when a cif file in CCP4 monomer library was given to phenix.refine. The cif file was $CLIBD_MON/3/3PG.cif. This was because phenix (mmtbx) ignores planarity restraint when dist_esd is 0. In the line 2055 of mmtbx/monomer_library/pdb_interpretation.py (phenix-1.10.1), plane_atom is not added to the list when (plane_atom.dist_esd in [None, 0]) is True. _chem_comp_plane_atom.dist_esd in 3PG.cif (ccp4-7.0) is all zero except the first atom in the plane definition. As far as I checked, there is no cif file where dist_esd of the first atom of plane is zero. At least, it would be helpful if phenix shows warning in this case. Best regards, Keitaro
Hi Keitaro,
My colleague noticed that planarity restraint of a ligand was ignored when a cif file in CCP4 monomer library was given to phenix.refine. The cif file was $CLIBD_MON/3/3PG.cif.
This was because phenix (mmtbx) ignores planarity restraint when dist_esd is 0. In the line 2055 of mmtbx/monomer_library/pdb_interpretation.py (phenix-1.10.1), plane_atom is not added to the list when (plane_atom.dist_esd in [None, 0]) is True.
yep, that's right. This is because restraints are T = SUM w*(p_model - p_ideal)**2, with w = 1/esd**2. So instead of bombing with division by zero error, the software ignores instances of restraints that would cause numeric troubles. I agree a warning may be helpful! Could you please send me a PDB file that contains this ligand? Pavel
Dear Pavel,
Yes, it would be problematic if 0 is not treated specially. What we
normally expect might be the planarity constraint if sigma is 0. I
don't know how ccp4 treats it (for example refmac5 keeps atoms planar
in this case), but maybe just the value of the first atom is used?
Thank you very much for your help but I am sorry that I cannot give a
PDB file because I do not have it.. Could you please generate a test
data yourself?
Best regards,
Keitaro
2016-09-28 0:23 GMT+09:00 Pavel Afonine
Hi Keitaro,
My colleague noticed that planarity restraint of a ligand was ignored when a cif file in CCP4 monomer library was given to phenix.refine. The cif file was $CLIBD_MON/3/3PG.cif.
This was because phenix (mmtbx) ignores planarity restraint when dist_esd is 0. In the line 2055 of mmtbx/monomer_library/pdb_interpretation.py (phenix-1.10.1), plane_atom is not added to the list when (plane_atom.dist_esd in [None, 0]) is True.
yep, that's right. This is because restraints are T = SUM w*(p_model - p_ideal)**2, with w = 1/esd**2. So instead of bombing with division by zero error, the software ignores instances of restraints that would cause numeric troubles.
I agree a warning may be helpful! Could you please send me a PDB file that contains this ligand?
Pavel
Hi Keitaro,
Yes, it would be problematic if 0 is not treated specially. What we normally expect might be the planarity constraint if sigma is 0.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea. What's done currently is there is a number, there is a formula that this number goes into, and there is a piece of code that uses the number and formula in a way to avoid numeric issues. As simple as this, no connotations involved. I believe using 0 as a trigger to tell that you want something else than what these numbers are used for otherwise is not the best way (why not put None or something?). All the best, Pavel
Dear Pavel, I have recently encountered a similar situation, where I took a CIF file for a ligand from the CCP4 monomer library and found planarity restraints are silently ignored in PHENIX. This is very confusing and dangerous.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea.
The "zero sigma" limit of a Gaussian is an infinitely sharp probability distribution (= delta function). So I think the transition from restraints to constraints is not "odd" but quite natural both mathematically and semantically. Anyway, if PHENIX does not want to treat 0 as constraints, I would prefer an error stop, not just a warning, because most people don't mean "ignore this record" by zero sigma. Takanori Nakane
It is my understanding that Refmac doesn’t make use of the CONST_X torsions in refinement either (I’d be very happy to hear more about this topic). In the monomer library files the planarities are defined by planes, so there is a redundancy in information. I think the CONST_X torsions are there for reporting but not for restraining (as they would be constraints by implication). If Phenix stopped because of torsions with 0 esd it would stop with all monomer library files with CONST_X records. We could report that they are there and not being used maybe. If someone defines an esd of 0 for a variable torsion then that might be a good place to throw an exception.
On Sep 27, 2016, at 6:19 PM, Takanori Nakane
wrote: Dear Pavel,
I have recently encountered a similar situation, where I took a CIF file for a ligand from the CCP4 monomer library and found planarity restraints are silently ignored in PHENIX. This is very confusing and dangerous.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea.
The "zero sigma" limit of a Gaussian is an infinitely sharp probability distribution (= delta function). So I think the transition from restraints to constraints is not "odd" but quite natural both mathematically and semantically.
Anyway, if PHENIX does not want to treat 0 as constraints, I would prefer an error stop, not just a warning, because most people don't mean "ignore this record" by zero sigma.
Takanori Nakane _______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb Unsubscribe: [email protected]
-- Paul Adams Division Director, Molecular Biophysics & Integrated Bioimaging, Lawrence Berkeley Lab Division Deputy for Biosciences, Advanced Light Source, Lawrence Berkeley Lab Adjunct Professor, Department of Bioengineering, U.C. Berkeley Vice President for Technology, the Joint BioEnergy Institute Laboratory Research Manager, ENIGMA Science Focus Area Building 33, Room 347 Building 80, Room 247 Building 978, Room 4126 Tel: 1-510-486-4225, Fax: 1-510-486-5909 http://cci.lbl.gov/paul Lawrence Berkeley Laboratory 1 Cyclotron Road BLDG 33R0345 Berkeley, CA 94720, USA. Executive Assistant: Louise Benvenue [ [email protected] ][ 1-510-495-2506 ] --
Hello Takanori, thanks for feedback! It is always very much appreciated!
I have recently encountered a similar situation, where I took a CIF file for a ligand from the CCP4 monomer library and found planarity restraints are silently ignored in PHENIX. This is very confusing and dangerous.
I can see how this may be confusing so I agree with Keitaro that a warning would be good to have. Also, phenix.refine always outputs .geo file that lists all the restraints for all atoms, which allows to check presence or absence of particular restraints for atoms in question.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea.
The "zero sigma" limit of a Gaussian is an infinitely sharp probability distribution (= delta function). So I think the transition from restraints to constraints is not "odd" but quite natural both mathematically and semantically.
Theoretically it is all right of course, but in terms of code it may be quite different pieces of code. It is much less cryptic if a decision to use one option (restraints) or another (constraints) is a made based on a clearly named parameter rather than inferring from interpreting a number.
Anyway, if PHENIX does not want to treat 0 as constraints, I would prefer an error stop, not just a warning,
It is always a balance between being conservative vs relaxed, automated vs less automated. There is no perfect solution for drawing the line, Iafraid. I'm sure if we make phenix.refine stop every time it sees esd=0 there will be a great deal of people flooding mailboxes saying "why it's stopping, I don't care, I want my refinement job finished", etc. phenix.refine favors automation paradigm while providing a lot of room for flexibility. You can define custom bonds, angles, torsions, planes, groups of planes to be parallel, as many as you wish and for any atoms of your choice. And you can make sure that what you defined was actually used by inspecting .geo files. A practical solution to "esd=0 issue" that is somewhere in between the two extremes ("stop if esd=0" vs "ignore and keep going") is to replace esd=0 with say esd=0.01 (or whatever is used for rings in TYR or PHE). This will let program going as before and you will get your plane that I bet one wouldn't be able to tell that it is not perfect just by looking at it!
because most people don't mean "ignore this record" by zero sigma.
Frankly, I don't know and I wouldn't be so sure. I'd say many wouldn't tell (and care about) the difference between restraints vs constrains, not to mention nuances about meaning of esd=0, which is fine. Those who don't mean "ignore this record" are probably trained to use the program X, where this assumption is made. Perhaps it's not too unreasonable to realize that not everything adopted in program X is equally and fully adopted in program Y. I can make a list of hundreds other similar decisions that are made under the hood of phenix.refine automatically. - Do hydrogens participate in bulk-solvent mask calculation? - Do vdw radii of non-H change when H are used in refinement? - What program does if it sees non-positive definite anisotropic ADP or any of T, L or S matrix? - What happens if input PDB has ANISOU records, TLS is not requested and resolution is not great for anisotropic refinement? - Are H treated as riding or freely refined if data resolution is 0.73A? - What data resolution is used when you request to do rigid-body refinement? - What refinement target is used (LS or ML) if there is less than 50 test (free) reflections? ... Yep, I can go on for 100+ more items like this. If every time we stop and ask users what to do that would not go well! All the best, Pavel
Hi Pavel, Thank you very much for clarification. I see your points and now agree that 'stop if esd is 0' is bit too much. Best regards, Takanori Nakane Am 2016年09月28日 um 11:33 schrieb Pavel Afonine:
Hello Takanori,
thanks for feedback! It is always very much appreciated!
I have recently encountered a similar situation, where I took a CIF file for a ligand from the CCP4 monomer library and found planarity restraints are silently ignored in PHENIX. This is very confusing and dangerous.
I can see how this may be confusing so I agree with Keitaro that a warning would be good to have. Also, phenix.refine always outputs .geo file that lists all the restraints for all atoms, which allows to check presence or absence of particular restraints for atoms in question.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea.
The "zero sigma" limit of a Gaussian is an infinitely sharp probability distribution (= delta function). So I think the transition from restraints to constraints is not "odd" but quite natural both mathematically and semantically.
Theoretically it is all right of course, but in terms of code it may be quite different pieces of code. It is much less cryptic if a decision to use one option (restraints) or another (constraints) is a made based on a clearly named parameter rather than inferring from interpreting a number.
Anyway, if PHENIX does not want to treat 0 as constraints, I would prefer an error stop, not just a warning,
It is always a balance between being conservative vs relaxed, automated vs less automated. There is no perfect solution for drawing the line, Iafraid. I'm sure if we make phenix.refine stop every time it sees esd=0 there will be a great deal of people flooding mailboxes saying "why it's stopping, I don't care, I want my refinement job finished", etc.
phenix.refine favors automation paradigm while providing a lot of room for flexibility. You can define custom bonds, angles, torsions, planes, groups of planes to be parallel, as many as you wish and for any atoms of your choice. And you can make sure that what you defined was actually used by inspecting .geo files.
A practical solution to "esd=0 issue" that is somewhere in between the two extremes ("stop if esd=0" vs "ignore and keep going") is to replace esd=0 with say esd=0.01 (or whatever is used for rings in TYR or PHE). This will let program going as before and you will get your plane that I bet one wouldn't be able to tell that it is not perfect just by looking at it!
because most people don't mean "ignore this record" by zero sigma.
Frankly, I don't know and I wouldn't be so sure. I'd say many wouldn't tell (and care about) the difference between restraints vs constrains, not to mention nuances about meaning of esd=0, which is fine.
Those who don't mean "ignore this record" are probably trained to use the program X, where this assumption is made. Perhaps it's not too unreasonable to realize that not everything adopted in program X is equally and fully adopted in program Y.
I can make a list of hundreds other similar decisions that are made under the hood of phenix.refine automatically. - Do hydrogens participate in bulk-solvent mask calculation? - Do vdw radii of non-H change when H are used in refinement? - What program does if it sees non-positive definite anisotropic ADP or any of T, L or S matrix? - What happens if input PDB has ANISOU records, TLS is not requested and resolution is not great for anisotropic refinement? - Are H treated as riding or freely refined if data resolution is 0.73A? - What data resolution is used when you request to do rigid-body refinement? - What refinement target is used (LS or ML) if there is less than 50 test (free) reflections? ... Yep, I can go on for 100+ more items like this. If every time we stop and ask users what to do that would not go well!
All the best, Pavel
_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb Unsubscribe: [email protected]
I've been following this from the sidelines with some interest! In terms of the point that a Gaussian with a zero sigma is a delta function, that is mathematically true but it would be asking a lot of programmers to take the limit of an integral in their evaluation code! Setting the sigma to an arbitrarily small number to simulate a constraint creates great potential numerical stability problems, because this makes both the first derivative (gradient) and second derivative (curvature) tend towards infinity. If you have an algorithm (like LBFGS) that infers the curvature from the behaviour of the gradient as parameters are varied, the algorithm typically can't cope with huge variations in curvature. Even if you're using a full second derivative matrix in Newton's method (which is rather uncommon in refinement programs), the matrix inversion step will become less numerically stable. So it seems that the best is to have sensible restraints for things that are allowed to vary and proper constraints for things that aren't allowed to vary. Whether or not you might choose to flag constraints by assigning zero sigmas is perhaps partly a matter of programming aesthetics, but it's much more obvious to the user that something very different will be done under the hood if the constraint flag is more explicit. Best wishes, Randy
On 28 Sep 2016, at 02:19, Takanori Nakane
wrote: Dear Pavel,
I have recently encountered a similar situation, where I took a CIF file for a ligand from the CCP4 monomer library and found planarity restraints are silently ignored in PHENIX. This is very confusing and dangerous.
with this logic you interpret the number and not use it. Follow this road you may end up with odd things, like let's treat 0 as "constraints", 0.1 as "strong restraints", 1.5 as "weak restraints", etc and in each case do something very different - I hope you get the idea.
The "zero sigma" limit of a Gaussian is an infinitely sharp probability distribution (= delta function). So I think the transition from restraints to constraints is not "odd" but quite natural both mathematically and semantically.
Anyway, if PHENIX does not want to treat 0 as constraints, I would prefer an error stop, not just a warning, because most people don't mean "ignore this record" by zero sigma.
Takanori Nakane _______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb Unsubscribe: [email protected]
------ Randy J. Read Department of Haematology, University of Cambridge Cambridge Institute for Medical Research Tel: + 44 1223 336500 Wellcome Trust/MRC Building Fax: + 44 1223 336827 Hills Road E-mail: [email protected] Cambridge CB2 0XY, U.K. www-structmed.cimr.cam.ac.uk
Here is a 3pg PDB file and a CCP4 restraints file. I should point out that Phenix applies the planar esd individually to each atom in the plane allowing some atoms to be more out of plane than others. I'm working on using this feature to better advantage at the moment. Cheers Nigel --- Nigel W. Moriarty Building 33R0349, Molecular Biophysics and Integrated Bioimaging Lawrence Berkeley National Laboratory Berkeley, CA 94720-8235 Phone : 510-486-5709 Email : [email protected] Fax : 510-486-5909 Web : CCI.LBL.gov On Tue, Sep 27, 2016 at 8:35 AM, Keitaro Yamashita < [email protected]> wrote:
Dear Pavel,
Yes, it would be problematic if 0 is not treated specially. What we normally expect might be the planarity constraint if sigma is 0. I don't know how ccp4 treats it (for example refmac5 keeps atoms planar in this case), but maybe just the value of the first atom is used?
Thank you very much for your help but I am sorry that I cannot give a PDB file because I do not have it.. Could you please generate a test data yourself?
Best regards, Keitaro
Hi Keitaro,
My colleague noticed that planarity restraint of a ligand was ignored when a cif file in CCP4 monomer library was given to phenix.refine. The cif file was $CLIBD_MON/3/3PG.cif.
This was because phenix (mmtbx) ignores planarity restraint when dist_esd is 0. In the line 2055 of mmtbx/monomer_library/pdb_interpretation.py (phenix-1.10.1), plane_atom is not added to the list when (plane_atom.dist_esd in [None, 0]) is True.
yep, that's right. This is because restraints are T = SUM w*(p_model - p_ideal)**2, with w = 1/esd**2. So instead of bombing with division by zero error, the software ignores instances of restraints that would cause numeric troubles.
I agree a warning may be helpful! Could you please send me a PDB file
2016-09-28 0:23 GMT+09:00 Pavel Afonine
: that contains this ligand?
Pavel
_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb Unsubscribe: [email protected]
participants (6)
-
Keitaro Yamashita
-
Nigel Moriarty
-
Paul Adams
-
Pavel Afonine
-
Randy Read
-
Takanori Nakane