From f00cee702083a13d69bcdf2d84cc47f2dfd1a127 Mon Sep 17 00:00:00 2001
From: "Nigel W. Moriarty" <nwmoriarty@lbl.gov>
Date: Wed, 16 Oct 2024 10:45:25 -0700
Subject: [PATCH] Working basic restraints validation

---
 mmtbx/monomer_library/geostd_utils.py         |  9 ++-
 .../programs/ligand_restraints_validation.py  | 68 +++++++++++++++++--
 2 files changed, 69 insertions(+), 8 deletions(-)

diff --git a/mmtbx/monomer_library/geostd_utils.py b/mmtbx/monomer_library/geostd_utils.py
index 7b193caf0e..b1c7235042 100644
--- a/mmtbx/monomer_library/geostd_utils.py
+++ b/mmtbx/monomer_library/geostd_utils.py
@@ -1,9 +1,9 @@
 from __future__ import division
 
+import iotbx
 from mmtbx.ligands.hierarchy_utils import _new_atom
 
 def get_as_hierarchy(filename):
-  import iotbx
   model = iotbx.cif.reader(filename).model()
   for key,block in model.items():
     if key=='comp_list':
@@ -26,9 +26,10 @@ def get_as_hierarchy(filename):
                        20.,
                        True,
                        )
+      # atom.set_serial(j+1)
       ag.append_atom(atom)
   rg = iotbx.pdb.hierarchy.residue_group()
-  rg.resseq='1'
+  rg.resseq='   1'
   rg.append_atom_group(ag)
   chain = iotbx.pdb.hierarchy.chain()
   chain.id='A'
@@ -39,3 +40,7 @@ def get_as_hierarchy(filename):
   ph.append_model(model)
   ph.reset_atom_i_seqs()
   return ph
+
+def as_cif_object(filename):
+  model = iotbx.cif.reader(filename).model()
+  return model
diff --git a/mmtbx/programs/ligand_restraints_validation.py b/mmtbx/programs/ligand_restraints_validation.py
index a930621e9f..3d2300c973 100644
--- a/mmtbx/programs/ligand_restraints_validation.py
+++ b/mmtbx/programs/ligand_restraints_validation.py
@@ -2,6 +2,8 @@
 
 from libtbx.program_template import ProgramTemplate
 
+# from cctbx.maptbx.box import shift_and_box_model
+
 class Program(ProgramTemplate):
 
   description = '''
@@ -27,15 +29,69 @@ def validate(self):
     pass
 
   # ---------------------------------------------------------------------------
+  def get_energies_sites(self, model, use_hydrogens=False):
+    assert model.restraints_manager is not None
+    if(use_hydrogens):
+      rm = model.restraints_manager
+      sc = model.get_sites_cart()
+    else:
+      hd_sel = model.get_xray_structure().hd_selection()
+      rm = model.restraints_manager.select(~hd_sel)
+      sc = model.get_sites_cart().select(~hd_sel)
+    energies_sites = \
+      rm.geometry.energies_sites(
+        sites_cart        = sc,
+        compute_gradients = False)
+    return energies_sites #.bond_deviations()[2]
+
   def processed_ligand_restaints(self, filename):
     from mmtbx.monomer_library.geostd_utils import get_as_hierarchy
+    from mmtbx.monomer_library.geostd_utils import as_cif_object
     hierarchy=get_as_hierarchy(filename)
-    hierarchy.show()
-    hierarchy.write_pdb_file('test.pdb')
-    print(hierarchy)
-    from mmtbx.model import model
-    m = model(pdb_hierarchy=hierarchy)
-    print(m)
+    hierarchy.sort_atoms_in_place()
+    # hierarchy.show()
+    starting_atoms=hierarchy.atoms()
+    starting_xyz=starting_atoms.extract_xyz()
+    if 0: hierarchy.write_pdb_file('test.pdb')
+    #
+    from mmtbx.model.model import manager
+    m = manager(pdb_hierarchy=hierarchy,
+                restraint_objects=[(filename, as_cif_object(filename))],
+                )
+    # m = shift_and_box_model(model = m, box_cushion = 5.)
+    # for atom in m.get_hierarchy().atoms(): print(atom.format_atom_record())
+    #
+    # geom min
+    #
+    from mmtbx.refinement import geometry_minimization
+    m.process(make_restraints=True)
+    geometry_minimization.run2(
+      restraints_manager = m.get_restraints_manager(),
+      pdb_hierarchy = m.get_hierarchy(),
+      correct_special_position_tolerance = 1.,
+      bond        = True,
+      nonbonded   = True,
+      angle       = True,
+      dihedral    = True,
+      chirality   = True,
+      planarity   = True,
+      parallelity = True,
+      )
+    # m.set_sites_cart_from_hierarchy()
+    # for atom in m.get_hierarchy().atoms(): print(atom.format_atom_record())
+    h=m.get_hierarchy()
+    if 0: h.write_pdb_file('test.pdb')
+    ending_xyz=h.atoms().extract_xyz()
+    self.results={}
+    rc = starting_xyz.rms_difference(ending_xyz)
+    print('\n  RMSD of starting and ending coordinates : %5.2f' % rc)
+    self.results['RMSD']=rc
+    es = self.get_energies_sites(m, use_hydrogens=True)
+    # es.show()
+    self.results['energies_sites']=es
+    rc = m.restraints_as_geo()
+    print(rc[:1000])
+    # print(rc)
 
   def run(self, log=None):
     for filename in self.data_manager.get_restraint_names():