geometric/internal.py
def remove_TR(self, xyz):
        """
        Project overall translation and rotation out of the DLCs.
        This feature is intended to be used when an optimization job appears
        to contain slow rotations of the whole system, which sometimes happens.
        Uses the same logic as build_dlc_1.
        """
        # Create three translation and three rotation primitive ICs for the whole system
        na = int(len(xyz)/3)
        alla = range(na)
        sel = xyz.reshape(-1,3).copy()
        TRPrims = []
        TRPrims.append(TranslationX(alla, w=np.ones(na)/na))
        TRPrims.append(TranslationY(alla, w=np.ones(na)/na))
        TRPrims.append(TranslationZ(alla, w=np.ones(na)/na))
        sel -= np.mean(sel, axis=0)
        rg = np.sqrt(np.mean(np.sum(sel**2, axis=1)))
        TRPrims.append(RotationA(alla, xyz, self.Prims.Rotators, w=rg))
        TRPrims.append(RotationB(alla, xyz, self.Prims.Rotators, w=rg))
        TRPrims.append(RotationC(alla, xyz, self.Prims.Rotators, w=rg))
        # If these primitive ICs are already there, then move them to the front
        primorder = []
        addPrims = []
        for prim in TRPrims:
            if prim in self.Prims.Internals:
                primorder.append(self.Prims.Internals.index(prim))
            else:
                addPrims.append(prim)
        for iprim, prim in enumerate(self.Prims.Internals):
            if prim not in TRPrims:
                primorder.append(iprim)
        self.Prims.Internals = addPrims + [self.Prims.Internals[p] for p in primorder]
        self.Vecs = np.vstack((np.zeros((len(addPrims), self.Vecs.shape[1]), dtype=float), self.Vecs[np.array(primorder), :]))
        
        self.clearCache()
        # Build DLCs with six extra in the front corresponding to the overall translations and rotations
        subVecs = self.Vecs.copy()
        self.Vecs = np.zeros((self.Vecs.shape[0], self.Vecs.shape[1]+6), dtype=float)
        self.Vecs[:6, :6] = np.eye(6)
        self.Vecs[:, 6:] = subVecs.copy()
        # pmat2d(self.Vecs, precision=3, format='f')
        # sys.exit()
        # This is the number of nonredundant DLCs that we expect to see
        G = self.Prims.GMatrix(xyz)
        Expect = np.sum(np.linalg.eigh(G)[0] > 1e-6)

        # Carry out Gram-Schmidt orthogonalization
        # Define a function for computing overlap
        def ov(vi, vj):
            return multi_dot([vi, G, vj])
        V = self.Vecs.copy()
        nv = V.shape[1]
        Vnorms = np.array([np.sqrt(ov(V[:, ic], V[:, ic])) for ic in range(nv)])
        # G_tmp = np.array([[ov(V[:, ic], V[:, jc]) for jc in range(nv)] for ic in range(nv)])
        # pmat2d(G_tmp, precision=3, format='f')
        # U holds the Gram-Schmidt orthogonalized DLCs
        U = np.zeros((V.shape[0], Expect), dtype=float)
        Unorms = np.zeros(Expect, dtype=float)
        ncon = len(self.Prims.cPrims)
        # Keep translations, rotations, and any constraints in sequential order
        # and project them out of the remaining DLCs
        for ic in range(6+ncon):
            U[:, ic] = V[:, ic].copy()
            ui = U[:, ic]
            Unorms[ic] = np.sqrt(ov(ui, ui))
            if Unorms[ic]/Vnorms[ic] < 0.1:
                logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original\n" % (ic-6, Unorms[ic]/Vnorms[ic]))
            V0 = V.copy()
            # Project out newest U column from all remaining V columns.
            for jc in range(ic+1, nv):
                vj = V[:, jc]
                vj -= ui * ov(ui, vj)/Unorms[ic]**2
        # Now keep the remaining DLC with the largest norm, perform projection,
        # then repeat until the expected number is found
        shift = 6+ncon
        for ic in range(shift, Expect):
            # Pick out the V column with the largest norm
            norms = np.array([np.sqrt(ov(V[:, jc], V[:, jc])) for jc in range(shift, nv)])
            imax = shift+np.argmax(norms)
            # Add this column to U
            U[:, ic] = V[:, imax].copy()
            ui = U[:, ic]
            Unorms[ic] = np.sqrt(ov(ui, ui))
            # Project out the newest U column from all V columns
            for jc in range(ncon, nv):
                V[:, jc] -= ui * ov(ui, V[:, jc])/Unorms[ic]**2
        # self.Vecs contains the linear combination coefficients that are our new DLCs
        self.Vecs = U[:, 6:].copy()
        self.Internals = ["Constraint" if i < ncon else "DLC" + " %i" % (i+1) for i in range(self.Vecs.shape[1])]

        Project overall translation and rotation out of the DLCs.
        This feature is intended to be used when an optimization job appears
        to contain slow rotations of the whole system, which sometimes happens.
        Uses the same logic as build_dlc_1.
        
Project overall translation and rotation out of the DLCs.
This feature is intended to be used when an optimization job appears
to contain slow rotations of the whole system, which sometimes happens.
Uses the same logic as build_dlc_1.
Project overall translation and rotation out of the DLCs . This feature is intended to be used when an optimization job appears to contain slow rotations of the whole system which sometimes happens . Uses the same logic as build_dlc_1 .
def remove_TR(self, xyz):
        na = int(len(xyz)/3)
        alla = range(na)
        sel = xyz.reshape(-1,3).copy()
        TRPrims = []
        TRPrims.append(TranslationX(alla, w=np.ones(na)/na))
        TRPrims.append(TranslationY(alla, w=np.ones(na)/na))
        TRPrims.append(TranslationZ(alla, w=np.ones(na)/na))
        sel -= np.mean(sel, axis=0)
        rg = np.sqrt(np.mean(np.sum(sel**2, axis=1)))
        TRPrims.append(RotationA(alla, xyz, self.Prims.Rotators, w=rg))
        TRPrims.append(RotationB(alla, xyz, self.Prims.Rotators, w=rg))
        TRPrims.append(RotationC(alla, xyz, self.Prims.Rotators, w=rg))
        primorder = []
        addPrims = []
        for prim in TRPrims:
            if prim in self.Prims.Internals:
                primorder.append(self.Prims.Internals.index(prim))
            else:
                addPrims.append(prim)
        for iprim, prim in enumerate(self.Prims.Internals):
            if prim not in TRPrims:
                primorder.append(iprim)
        self.Prims.Internals = addPrims + [self.Prims.Internals[p] for p in primorder]
        self.Vecs = np.vstack((np.zeros((len(addPrims), self.Vecs.shape[1]), dtype=float), self.Vecs[np.array(primorder), :]))
        self.clearCache()
        subVecs = self.Vecs.copy()
        self.Vecs = np.zeros((self.Vecs.shape[0], self.Vecs.shape[1]+6), dtype=float)
        self.Vecs[:6, :6] = np.eye(6)
        self.Vecs[:, 6:] = subVecs.copy()
        G = self.Prims.GMatrix(xyz)
        Expect = np.sum(np.linalg.eigh(G)[0] > 1e-6)
        def ov(vi, vj):
            return multi_dot([vi, G, vj])
        V = self.Vecs.copy()
        nv = V.shape[1]
        Vnorms = np.array([np.sqrt(ov(V[:, ic], V[:, ic])) for ic in range(nv)])
        U = np.zeros((V.shape[0], Expect), dtype=float)
        Unorms = np.zeros(Expect, dtype=float)
        ncon = len(self.Prims.cPrims)
        for ic in range(6+ncon):
            U[:, ic] = V[:, ic].copy()
            ui = U[:, ic]
            Unorms[ic] = np.sqrt(ov(ui, ui))
            if Unorms[ic]/Vnorms[ic] < 0.1:
                logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original\n" % (ic-6, Unorms[ic]/Vnorms[ic]))
            V0 = V.copy()
            for jc in range(ic+1, nv):
                vj = V[:, jc]
                vj -= ui * ov(ui, vj)/Unorms[ic]**2
        shift = 6+ncon
        for ic in range(shift, Expect):
            norms = np.array([np.sqrt(ov(V[:, jc], V[:, jc])) for jc in range(shift, nv)])
            imax = shift+np.argmax(norms)
            U[:, ic] = V[:, imax].copy()
            ui = U[:, ic]
            Unorms[ic] = np.sqrt(ov(ui, ui))
            for jc in range(ncon, nv):
                V[:, jc] -= ui * ov(ui, V[:, jc])/Unorms[ic]**2
        self.Vecs = U[:, 6:].copy()
        self.Internals = ["Constraint" if i < ncon else "DLC" + " %i" % (i+1) for i in range(self.Vecs.shape[1])]
def remove_TR ( self , xyz ) : na = int ( len ( xyz ) / 3 ) alla = range ( na ) sel = xyz . reshape ( - 1 , 3 ) . copy ( ) TRPrims = [ ] TRPrims . append ( TranslationX ( alla , w = np . ones ( na ) / na ) ) TRPrims . append ( TranslationY ( alla , w = np . ones ( na ) / na ) ) TRPrims . append ( TranslationZ ( alla , w = np . ones ( na ) / na ) ) sel -= np . mean ( sel , axis = 0 ) rg = np . sqrt ( np . mean ( np . sum ( sel ** 2 , axis = 1 ) ) ) TRPrims . append ( RotationA ( alla , xyz , self . Prims . Rotators , w = rg ) ) TRPrims . append ( RotationB ( alla , xyz , self . Prims . Rotators , w = rg ) ) TRPrims . append ( RotationC ( alla , xyz , self . Prims . Rotators , w = rg ) ) primorder = [ ] addPrims = [ ] for prim in TRPrims : if prim in self . Prims . Internals : primorder . append ( self . Prims . Internals . index ( prim ) ) else : addPrims . append ( prim ) for iprim , prim in enumerate ( self . Prims . Internals ) : if prim not in TRPrims : primorder . append ( iprim ) self . Prims . Internals = addPrims + [ self . Prims . Internals [ p ] for p in primorder ] self . Vecs = np . vstack ( ( np . zeros ( ( len ( addPrims ) , self . Vecs . shape [ 1 ] ) , dtype = float ) , self . Vecs [ np . array ( primorder ) , : ] ) ) self . clearCache ( ) subVecs = self . Vecs . copy ( ) self . Vecs = np . zeros ( ( self . Vecs . shape [ 0 ] , self . Vecs . shape [ 1 ] + 6 ) , dtype = float ) self . Vecs [ : 6 , : 6 ] = np . eye ( 6 ) self . Vecs [ : , 6 : ] = subVecs . copy ( ) G = self . Prims . GMatrix ( xyz ) Expect = np . sum ( np . linalg . eigh ( G ) [ 0 ] > 1e-6 ) def ov ( vi , vj ) : return multi_dot ( [ vi , G , vj ] ) V = self . Vecs . copy ( ) nv = V . shape [ 1 ] Vnorms = np . array ( [ np . sqrt ( ov ( V [ : , ic ] , V [ : , ic ] ) ) for ic in range ( nv ) ] ) U = np . zeros ( ( V . shape [ 0 ] , Expect ) , dtype = float ) Unorms = np . zeros ( Expect , dtype = float ) ncon = len ( self . Prims . cPrims ) for ic in range ( 6 + ncon ) : U [ : , ic ] = V [ : , ic ] . copy ( ) ui = U [ : , ic ] Unorms [ ic ] = np . sqrt ( ov ( ui , ui ) ) if Unorms [ ic ] / Vnorms [ ic ] < 0.1 : logger . warning ( "Constraint %i is almost redundant; after projection norm is %.3f of original\n" % ( ic - 6 , Unorms [ ic ] / Vnorms [ ic ] ) ) V0 = V . copy ( ) for jc in range ( ic + 1 , nv ) : vj = V [ : , jc ] vj -= ui * ov ( ui , vj ) / Unorms [ ic ] ** 2 shift = 6 + ncon for ic in range ( shift , Expect ) : norms = np . array ( [ np . sqrt ( ov ( V [ : , jc ] , V [ : , jc ] ) ) for jc in range ( shift , nv ) ] ) imax = shift + np . argmax ( norms ) U [ : , ic ] = V [ : , imax ] . copy ( ) ui = U [ : , ic ] Unorms [ ic ] = np . sqrt ( ov ( ui , ui ) ) for jc in range ( ncon , nv ) : V [ : , jc ] -= ui * ov ( ui , V [ : , jc ] ) / Unorms [ ic ] ** 2 self . Vecs = U [ : , 6 : ] . copy ( ) self . Internals = [ "Constraint" if i < ncon else "DLC" + " %i" % ( i + 1 ) for i in range ( self . Vecs . shape [ 1 ] ) ]
Project overall translation and rotation out of the DLCs.
Project overall translation and rotation out of the DLCs .
"""
        Project overall translation and rotation out of the DLCs.
        This feature is intended to be used when an optimization job appears
        to contain slow rotations of the whole system, which sometimes happens.
        Uses the same logic as build_dlc_1.
        """ # Create three translation and three rotation primitive ICs for the whole system # If these primitive ICs are already there, then move them to the front # Build DLCs with six extra in the front corresponding to the overall translations and rotations # pmat2d(self.Vecs, precision=3, format='f') # sys.exit() # This is the number of nonredundant DLCs that we expect to see # Carry out Gram-Schmidt orthogonalization # Define a function for computing overlap # G_tmp = np.array([[ov(V[:, ic], V[:, jc]) for jc in range(nv)] for ic in range(nv)]) # pmat2d(G_tmp, precision=3, format='f') # U holds the Gram-Schmidt orthogonalized DLCs # Keep translations, rotations, and any constraints in sequential order # and project them out of the remaining DLCs # Project out newest U column from all remaining V columns. # Now keep the remaining DLC with the largest norm, perform projection, # then repeat until the expected number is found # Pick out the V column with the largest norm # Add this column to U # Project out the newest U column from all V columns # self.Vecs contains the linear combination coefficients that are our new DLCs
self
xyz