diff --git a/mmtbx/geometry_restraints/mopac_manager.py b/mmtbx/geometry_restraints/mopac_manager.py index b904603609..22add0a6f3 100644 --- a/mmtbx/geometry_restraints/mopac_manager.py +++ b/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) diff --git a/mmtbx/geometry_restraints/quantum_interface.py b/mmtbx/geometry_restraints/quantum_interface.py index 008d4aef8d..a517d22959 100644 --- a/mmtbx/geometry_restraints/quantum_interface.py +++ b/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): diff --git a/mmtbx/geometry_restraints/quantum_restraints_manager.py b/mmtbx/geometry_restraints/quantum_restraints_manager.py index 2464474791..27b3476c5a 100644 --- a/mmtbx/geometry_restraints/quantum_restraints_manager.py +++ b/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,9 +534,11 @@ 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, @@ -535,6 +546,13 @@ def get_qm_manager(ligand_model, buffer_model, qmr, program_goal, log=StringIO() 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, diff --git a/mmtbx/ligands/electrons.py b/mmtbx/ligands/electrons.py index efc417c008..72b3e91c9f 100644 --- a/mmtbx/ligands/electrons.py +++ b/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)] diff --git a/mmtbx/model/model.py b/mmtbx/model/model.py index 1db3a0007a..a1db5d5ffd 100644 --- a/mmtbx/model/model.py +++ b/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()