better QM multiplicity support

mmtbx/geometry_restraints/mopac_manager.py

@@ -6,6 +6,7 @@
 from mmtbx.geometry_restraints import base_qm_manager
 
 import libtbx.load_env
+from libtbx import Auto
 
 from mmtbx.geometry_restraints.base_qm_manager import get_internal_coordinate_value
 
@@ -34,11 +35,22 @@ def _input_header(self):
       nproc_str=''
     else:
       nproc_str='THREADS=%s' % self.nproc
-    outl = '%s %s %s %s DISP \n%s\n\n' % (
+    multiplicity_str=''
+    if self.multiplicity not in [None, Auto]:
+      multiplicity_str=[None,
+                        'singlet', # - 0 unpaired electrons
+                        'doublet', # - 1 unpaired electrons
+                        'triplet', # - 2 unpaired electrons
+                        'quartet', # - 3 unpaired electrons
+                        'quintet', # - 4 unpaired electrons
+                        'sextet', # - 5 unpaired electrons
+                        ][self.multiplicity]
+    outl = '%s %s %s %s DISP %s\n%s\n\n' % (
      self.method,
      self.basis_set,
      self.solvent_model,
      'CHARGE=%s %s' % (self.charge, nproc_str),
+     multiplicity_str,
      self.preamble,
      )
     return outl
@@ -120,11 +132,12 @@ def get_input_lines(self, optimise_ligand=True, optimise_h=True, constrain_torsi
   def get_input_filename(self):
     return 'mopac_%s.mop' % self.preamble
 
-  def read_xyz_output(self):
+  def read_xyz_output(self, verbose=False):
     filename = self.get_coordinate_filename()
     filename = filename.replace('.arc', '.out')
     if not os.path.exists(filename):
       raise Sorry('QM output filename not found: %s' % filename)
+    if verbose: print('reading %s' % filename)
     f=open(filename, 'r')
     lines = f.read()
     del f
@@ -179,6 +192,36 @@ def read_charge(self, filename=None):
         break
     return self.charge
 
+  def validate_atomic_charges(self, filename=None):
+    '''
+                 NET ATOMIC CHARGES AND DIPOLE CONTRIBUTIONS
+
+  ATOM NO.   TYPE          CHARGE      No. of ELECS.   s-Pop       p-Pop
+    1          C          -0.691089        4.6911     1.10997     3.58112
+    2          C           0.963519        3.0365     1.15413     1.88235
+    3          O          -0.111514        6.1115     1.65968     4.45183
+    4          H           0.405761        0.5942     0.59424
+    5          H           0.399617        0.6004     0.60038
+    6          H           0.404684        0.5953     0.59532
+    7          H           1.000000        0.0000     0.00000
+    8          H           1.000000        0.0000     0.00000
+    9          H           0.629022        0.3710     0.37098
+ DIPOLE           X         Y         Z       TOTAL
+ POINT-CHG.    29.803  -147.398  -116.121   189.996
+ HYBRID         0.786     0.431     0.475     1.015
+ SUM           30.590  -146.966  -115.646   189.496'''
+    lines = self.get_lines(filename=filename)
+    reading=False
+    for line in lines.splitlines():
+      if reading:
+        print(line)
+        assert 0
+      if line.find('ATOM NO.   TYPE          CHARGE      No. of ELECS.')>-1:
+        reading=True
+
+  def validate_calculation(self):
+    self.validate_atomic_charges()
+
   def read_gradients(self, filename=None):
     lines = self.get_lines(filename=filename)
     print(lines)

mmtbx/geometry_restraints/quantum_interface.py

@@ -84,7 +84,20 @@ def get_qm_restraints_scope(validate=True, verbose=False):
       .type = atom_selection
       .input_size = 400
     charge = None
-      .type = float
+      .type = int
+  }
+  specific_atom_multiplicities
+    .optional = True
+    .multiple = True
+    .short_caption = Specify the multiplicity for a specific atom (mostly metal ions). \
+    Selection is not needed but is great for bookkeeping.
+    .style = auto_align
+  {
+    atom_selection = None
+      .type = atom_selection
+      .input_size = 400
+    multiplicity = None
+      .type = int
   }
 
   calculate = *in_situ_opt starting_energy final_energy \
@@ -183,13 +196,17 @@ def get_qm_restraints_scope(validate=True, verbose=False):
 
 master_phil_str = get_qm_restraints_scope()
 
-def electrons(model, specific_atom_charges=None, log=None):
+def electrons(model,
+              specific_atom_charges=None,
+              specific_atom_multiplicities=None,
+              log=None):
   from libtbx.utils import Sorry
   from mmtbx.ligands import electrons
   atom_valences = electrons.electron_distribution(
     model.get_hierarchy(), # needs to be altloc free
     model.get_restraints_manager().geometry,
     specific_atom_charges=specific_atom_charges,
+    specific_atom_multiplicities=specific_atom_multiplicities,
     log=log,
     verbose=False,
   )
@@ -208,9 +225,14 @@ def electrons(model, specific_atom_charges=None, log=None):
         print('    %s' % i[0], file=log)
     # raise Sorry('Unusual charges found')
   charged_atoms = atom_valences.get_charged_atoms()
+  print('''
+  Complete valence picture
+%s
+  '''% (atom_valences.show()), file=log)
   return atom_valences.get_total_charge()
 
 def get_safe_filename(s, compact_selection_syntax=True):
+  import string
   assert compact_selection_syntax
   if compact_selection_syntax:
     s=s.replace('chain', '')
@@ -222,12 +244,18 @@ def get_safe_filename(s, compact_selection_syntax=True):
     s=s.replace('  ',' ')
   if s[0]==' ': s=s[1:]
   s=s.replace(' ','_')
+  for i in range(26):
+    a=string.ascii_uppercase[i]
+    s=s.replace("'%s'"%a, a)
   s=s.replace("'",'_prime_')
   s=s.replace('*','_star_')
   s=s.replace('(','_lb_')
   s=s.replace(')','_rb_')
   s=s.replace('=', '_equals_')
   s=s.replace(':', '_colon_')
+  s=s.replace('"', '_quote_')
+  while s.find('__')>-1:
+    s=s.replace('__','_')
   return s
 
 def populate_qmr_defaults(qmr):
@@ -266,6 +294,13 @@ def get_working_directory(model, params, prefix=None):
     rc='%s_%s' % (prefix, rc)
   return rc
 
+def get_total_multiplicity(qmr):
+  tm = 0
+  macs = qmr.specific_atom_multiplicities
+  for mac in macs:
+    tm += mac.multiplicity
+  return max(tm,1)
+
 def get_preamble(macro_cycle, i, qmr, old_style=False, compact_selection_syntax=True):
   qmr = populate_qmr_defaults(qmr)
   s=''
@@ -297,6 +332,9 @@ def get_preamble(macro_cycle, i, qmr, old_style=False, compact_selection_syntax=
     s+='_%s' % get_safe_filename(qmr.package.basis_set)
   if qmr.package.solvent_model is not Auto and qmr.package.solvent_model:
     s+='_%s' % get_safe_filename(qmr.package.solvent_model)
+  multiplicity=get_total_multiplicity(qmr)
+  if multiplicity!=1:
+    s+='_%s' % (multiplicity)
   return s
 
 def is_any_quantum_package_installed(env):

mmtbx/geometry_restraints/quantum_restraints_manager.py

@@ -185,7 +185,7 @@ def _reindexing(mod, sel, verbose=False):
       atoms[ra[i_seq]].tmp=i_seq
     if verbose:
       for atom in mod.get_atoms():
-        print(atom.quote(), atom.i_seq, atom.tmp)
+        print('...',atom.quote(), atom.i_seq, atom.tmp)
     return mod
   assert (selection_array, selection_str).count(None)==1
   if selection_str:
@@ -394,8 +394,17 @@ def get_ligand_buffer_models(model, qmr, verbose=False, write_steps=False, log=N
   if debug: print('buffer_selection_string',buffer_selection_string)
   buffer_model = select_and_reindex(model, buffer_selection_string)
   if write_steps: write_pdb_file(buffer_model, 'pre_remove_altloc.pdb', None)
-  if 0: retain_only_one_alternative_conformation(buffer_model, 'B')
-  buffer_model.remove_alternative_conformations(always_keep_one_conformer=True)
+  altloc=None
+  if qmr.selection.find('altloc')>-1:
+    tmp = qmr.selection.split()
+    i = tmp.index('altloc')
+    altloc=tmp[i+1]
+    altloc=altloc.replace(')','')
+    altloc=altloc.replace('"','')
+    altloc=altloc.replace("'",'')
+  if 0: retain_only_one_alternative_conformation(buffer_model, altloc)
+  buffer_model.remove_alternative_conformations(always_keep_one_conformer=True,
+                                                altloc_to_keep=altloc)
   if write_steps: write_pdb_file(buffer_model, 'post_remove_altloc.pdb', None)
   validate_ligand_buffer_models(ligand_model, buffer_model, qmr, log=log)
   if write_steps: write_pdb_file(buffer_model, 'pre_super_cell.pdb', None)
@@ -525,16 +534,25 @@ def get_qm_manager(ligand_model, buffer_model, qmr, program_goal, log=StringIO()
   else:
     assert 0, 'program_goal %s not in list' % program_goal
   specific_atom_charges = qmr.specific_atom_charges
+  specific_atom_multiplicities = qmr.specific_atom_multiplicities
   total_charge = quantum_interface.electrons(
     electron_model,
     specific_atom_charges=specific_atom_charges,
+    specific_atom_multiplicities=specific_atom_multiplicities,
     log=log)
   if total_charge!=qmr.package.charge:
     print(u'  Update charge for "%s" cluster : %s ~> %s' % (qmr.selection,
                                                             qmr.package.charge,
                                                             total_charge),
           file=log)
     qmr.package.charge=total_charge
+  total_multiplicity = quantum_interface.get_total_multiplicity(qmr)
+  if total_multiplicity!=qmr.package.multiplicity:
+    print(u'  Update multiplicity for "%s" cluster : %s ~> %s' % (qmr.selection,
+                                                                  qmr.package.multiplicity,
+                                                                  total_multiplicity),
+          file=log)
+    qmr.package.multiplicity=total_multiplicity
   qmm = qmm(electron_model.get_atoms(),
             qmr.package.method,
             qmr.package.basis_set,

mmtbx/ligands/electrons.py

@@ -105,7 +105,7 @@ def __init__(self,
                hierarchy,
                grm,
                specific_atom_charges=None, # a list of selections and charges
-               specific_atom_multicities=None,
+               specific_atom_multiplicities=None,
                alternative_location_id=None,
                alternative_location_index=None,
                log=None,
@@ -123,6 +123,8 @@ def __init__(self,
       sites_cart=self.xrs.sites_cart())
     #
     self.specific_atom_charges = specific_atom_charges
+    self.specific_atom_multiplicities = specific_atom_multiplicities
+    self.atoms_with_charges_set = []
     if log is None:
       self.logger = sys.stdout
     else:
@@ -147,6 +149,7 @@ def __init__(self,
     self.set_charges() # place for selections
     self.validate_metals()
     self.form_bonds()
+    self.adjust_for_multiplicity()
 
   def __repr__(self):
     return self.show()
@@ -278,8 +281,24 @@ def set_charges(self):
         metal_asc = self.hierarchy.atom_selection_cache()
         metal_sel = metal_asc.selection(sac.atom_selection)
         metal_hierarchy = self.hierarchy.select(metal_sel)
-        for atom in metal_hierarchy.atoms():
+        for i, atom in enumerate(metal_hierarchy.atoms()):
           self[atom.i_seq]=sac.charge*-1
+          self.atoms_with_charges_set.append(atom.i_seq)
+          assert i<1
+
+  def adjust_for_multiplicity(self):
+    if self.specific_atom_multiplicities:
+      for i, mac in enumerate(self.specific_atom_multiplicities):
+        radical_asc = self.hierarchy.atom_selection_cache()
+        radical_sel = radical_asc.selection(mac.atom_selection)
+        radical_hierarchy = self.hierarchy.select(radical_sel)
+        for i, atom in enumerate(radical_hierarchy.atoms()):
+          if self[atom.i_seq] and not atom.i_seq in self.atoms_with_charges_set:
+            if mac.multiplicity==2:
+              self[atom.i_seq]=0
+              print('\nCharge on %s changed to zero because of multiplicity. CHECK!' % atom.quote(),
+                    file=self.logger)
+              break
 
   def _has_metal(self, atom1, atom2):
     is_metal_count = [0,1][self.properties.is_metal(atom1.element)]

mmtbx/model/model.py

@@ -3669,13 +3669,14 @@ def show_groups(self, rigid_body = None, tls = None,
     out.flush()
     time_model_show += timer.elapsed()
 
-  def remove_alternative_conformations(self, always_keep_one_conformer):
+  def remove_alternative_conformations(self, always_keep_one_conformer, altloc_to_keep=None):
     # XXX This is not working correctly when something was deleted.
     # Need to figure out a way to update everything so GRM
     # construction will not fail.
     self.geometry_restraints = None
     self._pdb_hierarchy.remove_alt_confs(
-        always_keep_one_conformer=always_keep_one_conformer)
+        always_keep_one_conformer=always_keep_one_conformer,
+        altloc_to_keep=altloc_to_keep)
     self._pdb_hierarchy.sort_atoms_in_place()
     self._pdb_hierarchy.atoms_reset_serial()
     self.update_xrs()