Merge pull request #120 from molssi-seamm/dev

Supporting space group symmetry

HISTORY.rst

@@ -1,6 +1,8 @@
 =======
 History
 =======
+2023.8.30 -- Support for spacegroup symmetry
+
 2023.7.27 -- Bugfix: printing bond order info
  * If the bond orders were printed but not used on the system, the code crashed.
 

mopac_step/energy.py

@@ -8,6 +8,7 @@
 from pathlib import Path
 import textwrap
 
+import numpy as np
 from tabulate import tabulate
 
 import mopac_step
@@ -903,15 +904,26 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
  directory.mkdir(parents=True, exist_ok=True)
 
  system, configuration = self.get_system_configuration(None)
- symbols = configuration.atoms.symbols
+ symbols = configuration.atoms.asymmetric_symbols
  atoms = configuration.atoms
+ symmetry = configuration.symmetry
  if "ATOM_CHARGES" in data:
  # Add to atoms (in coordinate table)
  if "charge" not in atoms:
  atoms.add_attribute(
  "charge", coltype="float", configuration_dependent=True
  )
- atoms["charge"][0:] = data["ATOM_CHARGES"]
+ if symmetry.n_symops == 1:
+ chgs = data["ATOM_CHARGES"]
+ else:
+ chgs, delta = symmetry.symmetrize_atomic_scalar(data["ATOM_CHARGES"])
+ delta = np.array(delta)
+ max_delta = np.max(abs(delta))
+ text_lines.append(
+ "The maximum difference of the charges of symmetry related atoms "
+ f"was {max_delta:.4f}\n"
+ )
+ atoms["charge"][0:] = chgs
 
  # Print the charges and dump to a csv file
  chg_tbl = {
@@ -932,15 +944,25 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
  atoms.add_attribute(
  "spin", coltype="float", configuration_dependent=True
  )
- atoms["spin"][0:] = spins
+ if symmetry.n_symops == 1:
+ atoms["spin"][0:] = spins
+ else:
+ spins, delta = symmetry.symmetrize_atomic_scalar(spins)
+ atoms["spins"][0:] = spins
+ delta = np.array(delta)
+ max_delta = np.max(abs(delta))
+ text_lines.append(
+ " The maximum difference of the spins of symmetry "
+ f"related atoms was {max_delta:.4f}.\n"
+ )
 
  header = " Atomic charges and spins"
  chg_tbl["Spin"] = []
  writer.writerow(["Atom", "Element", "Charge", "Spin"])
  for atom, symbol, q, s in zip(
  range(1, len(symbols) + 1),
  symbols,
- data["ATOM_CHARGES"],
+ chgs,
  spins,
  ):
  q = f"{q:.3f}"
@@ -956,7 +978,7 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
  for atom, symbol, q in zip(
  range(1, len(symbols) + 1),
  symbols,
- data["ATOM_CHARGES"],
+ chgs,
  ):
  q = f"{q:.2f}"
  writer.writerow([atom, symbol, q])
@@ -1047,17 +1069,7 @@ def _bond_orders(self, control, bond_order_matrix, configuration):
  options = self.parent.options
  text_lines = []
  if len(symbols) <= int(options["max_atoms_to_print"]):
- if "name" in configuration.atoms:
- name = configuration.atoms.get_column_data("name")
- else:
- name = []
- count = {}
- for symbol in symbols:
- if symbol not in count:
- count[symbol] = 1
- else:
- count[symbol] += 1
- name.append(f"{symbol}{count[symbol]}")
+ name = configuration.atoms.names
  table = {
  "i": [name[i] for i in bond_i],
  "j": [name[j] for j in bond_j],

mopac_step/optimization.py

@@ -7,6 +7,8 @@
 import textwrap
 import traceback
 
+import numpy as np
+
 import mopac_step
 import seamm
 import seamm_util.printing as printing
@@ -67,14 +69,14 @@ def description_text(self, P=None):
  text += " The cell for periodic systems will not be optimized."
 
  if P["gnorm"] == "default":
- text += " The geometrical convergence is the default of " "1.0 kcal/mol/Å."
+ text += " The geometrical convergence is the default of 1.0 kcal/mol/Å."
  elif P["gnorm"][0] == "$":
  text += (
  " The geometrical convergence will be determined "
  "at runtime by '{gnorm}'."
  )
  else:
- text += " The geometrical convergence is {gnorm} kcal/mol/Å."
+ text += " The geometrical convergence is {gnorm}."
 
  # Put in the description of the energy calculation
  text += "\n\nThe energy and forces will be c" + energy_description[1:]
@@ -399,18 +401,24 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
  periodicity=periodicity,
  atomset=starting_configuration.atomset,
  cell_id=starting_configuration.cell_id,
+ symmetry=starting_configuration.symmetry_id,
  )
  else:
  configuration = system.create_configuration(
  periodicity=periodicity,
  atomset=starting_configuration.atomset,
  bondset=starting_configuration.bondset,
  cell_id=starting_configuration.cell_id,
+ symmetry=starting_configuration.symmetry_id,
  )
  configuration.charge = starting_configuration.charge
  configuration.spin_multiplicity = (
  starting_configuration.spin_multiplicity
  )
+ # Add initial coordinates so symmetry is handled properly
+ configuration.atoms.set_coordinates(
+ starting_configuration.atoms.get_coordinates(asymmetric=True)
+ )
  else:
  configuration = starting_configuration
 
@@ -442,6 +450,19 @@ def analyze(self, indent="", data_sections=[], out_sections=[], table=None):
  it = iter(data["ATOM_X_UPDATED"][-1])
  for x in it:
  xyz.append([float(x), float(next(it)), float(next(it))])
+
+ if configuration.symmetry.n_symops > 1:
+ # Convert to coordinates of just the asymmetric atoms.
+ xyz, delta = configuration.symmetry.symmetrize_coordinates(
+ xyz, fractionals=False
+ )
+ delta = np.array(delta)
+ displacement = np.linalg.norm(delta, axis=1)
+ max_disp = displacement.max()
+ text += (
+ f"\nThe largest displacement from symmetry was {max_disp:.6f} Å.\n"
+ )
+
  configuration.atoms.set_coordinates(xyz, fractionals=False)
 
  # And the name of the configuration.