Skip to content

Commit

Permalink
Ensure rigid optimization works when one of the fragments is monatomi…
Browse files Browse the repository at this point in the history
…c or diatomic
  • Loading branch information
leeping committed Apr 11, 2024
1 parent 7789a53 commit 6a152e8
Show file tree
Hide file tree
Showing 2 changed files with 9 additions and 7 deletions.
14 changes: 8 additions & 6 deletions geometric/internal.py
Original file line number Diff line number Diff line change
Expand Up @@ -302,7 +302,8 @@ def second_derivative(self, xyz):

def IsTR(Prim):
return type(Prim) in (TranslationX, TranslationY, TranslationZ,
RotationA, RotationB, RotationC)
RotationA, RotationB, RotationC,
CartesianX, CartesianY, CartesianZ)

class TranslationX(PrimitiveCoordinate):
def __init__(self, a, w):
Expand Down Expand Up @@ -2167,13 +2168,14 @@ def rigid(self, val):
self._rigid = val

class PrimitiveInternalCoordinates(InternalCoordinates):
def __init__(self, molecule, connect=False, addcart=False, constraints=None, cvals=None, **kwargs):
def __init__(self, molecule, connect=False, addcart=False, constraints=None, cvals=None, connect_isolated=True, **kwargs):
super(PrimitiveInternalCoordinates, self).__init__()
# connect = True corresponds to "traditional" internal coordinates with minimum spanning bonds
self.connect = connect
# connect = False, addcart = True corresponds to HDLC
# connect = False, addcart = False corresponds to TRIC
self.addcart = addcart
self.connect_isolated = connect_isolated
if 'rigid' in kwargs and kwargs['rigid'] is not None:
raise RuntimeError('Do not use rigid molecules with PrimitiveInternalCoordinates')
self.Internals = []
Expand All @@ -2197,7 +2199,6 @@ def makePrimitives(self, molecule, connect, addcart):
# force_bonds=False is set because we don't want to override
# bond order-based bonds that may have been obtained earlier.
molecule.build_topology(force_bonds=False)
connect_isolated = True
if 'resid' in molecule.Data.keys():
# Create fragments corresponding to unique resID numbers if provided.
frags = []
Expand All @@ -2218,7 +2219,7 @@ def makePrimitives(self, molecule, connect, addcart):
else:
# Create fragments based on connectivity in provided molecule object.
frags = [list(m.nodes()) for m in molecule.molecules]
if connect_isolated:
if self.connect_isolated:
isolates = [list(m.nodes())[0] for m in molecule.molecules if len(m.nodes()) == 1]
for i in isolates:
j = molecule.get_closest_atom(i, pbc=False)[0]
Expand Down Expand Up @@ -2963,7 +2964,7 @@ def covalent(a, b):


class DelocalizedInternalCoordinates(InternalCoordinates):
def __init__(self, molecule, imagenr=0, build=False, connect=False, addcart=False, constraints=None, cvals=None, rigid=False, remove_tr=False, cart_only=False, conmethod=0):
def __init__(self, molecule, imagenr=0, build=False, connect=False, addcart=False, constraints=None, cvals=None, rigid=False, remove_tr=False, cart_only=False, conmethod=0, connected_isolated=True):
super(DelocalizedInternalCoordinates, self).__init__()
# cart_only is just because of how I set up the class structure.
if cart_only: return
Expand All @@ -2984,7 +2985,8 @@ def __init__(self, molecule, imagenr=0, build=False, connect=False, addcart=Fals
if connect or addcart:
raise RuntimeError("Rigid optimizations not available using non-TRIC coordinates")
# The DLC contains an instance of primitive internal coordinates.
self.Prims = PrimitiveInternalCoordinates(molecule, connect=connect, addcart=addcart, constraints=constraints, cvals=cvals)
if rigid: connect_isolated = False
self.Prims = PrimitiveInternalCoordinates(molecule, connect=connect, addcart=addcart, constraints=constraints, cvals=cvals, connect_isolated=connect_isolated)
self.frags = self.Prims.frags
self.na = molecule.na
# Atomic mass array
Expand Down
2 changes: 1 addition & 1 deletion geometric/molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -1906,7 +1906,7 @@ def atom_select(self,atomslice,build_topology=True):
New.Data['bonds'] = [(list(atomslice).index(b[0]), list(atomslice).index(b[1])) for b in self.bonds if (b[0] in atomslice and b[1] in atomslice)]
New.top_settings = copy.deepcopy(self.top_settings)

if build_topology:
if New.na > 1 and build_topology:
New.build_topology(force_bonds=False)
return New

Expand Down

0 comments on commit 6a152e8

Please sign in to comment.