diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 559ff36..45760d5 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -32,6 +32,7 @@ jobs: shell: bash -l {0} run: | conda install -y pytest pytest-cov + conda install -y pip - name: Install MolSym shell: bash -l {0} run: | diff --git a/README.md b/README.md index 7d83bd0..e4c9731 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# ConcordantModes +# concordantmodes

Code style: black @@ -17,7 +17,7 @@ The following procedure may be used to enact the CMA-0A procedure using manual s ## Installation User must have at least: -Python version 3.7 +Python version 3.10 Numpy version 1.13.3 qcelemental @@ -27,7 +27,7 @@ Pytest-xdist 3.0.2 (-c condaforge) A simple way to ensure this program may be run is by creating a conda environment with the following command: -`conda create --name CMA python=3.9 sympy numpy scipy qcelemental` +`conda create --name CMA python=3.10 sympy numpy scipy qcelemental` If you need to install anaconda, consult the [official website](https://www.anaconda.com/products/individual) diff --git a/concordantmodes/algorithm.py b/concordantmodes/algorithm.py index 4059117..f87a521 100644 --- a/concordantmodes/algorithm.py +++ b/concordantmodes/algorithm.py @@ -4,11 +4,85 @@ from scipy import stats -class Algorithm(object): +class Algorithm: """ - The purpose of this class is to return a list of indices by which the force constants of the CMA method - will be computed. These indices will be determined by user input where nonabelian symmetry can be - exploited. + Generate force-constant displacement indices for the Concordant Modes + Algorithm (CMA). + + This class determines which coordinate pairs `(i, j)` will be evaluated + during finite-difference force-constant calculations. The generated index + list depends on: + + - The CMA level being performed (Level A or Level B). + - Whether force constants are computed from energies or gradients. + - Whether molecular symmetry is enabled. + - The irreducible-representation (irrep) structure of the projected + vibrational coordinate space. + + For Level B calculations without gradient-based Hessian construction, + off-diagonal force constants may be included. For Level A calculations + (or when gradient derivatives are used), only diagonal force constants + are generated. When symmetry is enabled, indices are restricted to + coordinate pairs belonging to the same irreducible representation, + reducing the number of required electronic structure calculations. + + Parameters + ---------- + num_deg_free : int + Number of vibrational degrees of freedom in the projected coordinate + space. + cma_level : {"A", "B"} + CMA stage for which displacement indices are generated. + + - ``"B"``: Lower-level reference Hessian calculation. + - ``"A"``: Higher-level CMA correction calculation. + + options : object + User options containing symmetry settings and derivative-level + controls. The following attributes are used: + + - ``molsym_symmetry`` : bool + - ``deriv_level_b`` : bool + + proj_irreps : list, optional + Symmetry decomposition of the projected coordinate space. Each + element specifies the number of coordinates belonging to an + irreducible representation. Degenerate irreducible + representations may be represented as nested lists. + + Attributes + ---------- + indices : list[list[int]] + Flat list of coordinate index pairs to be displaced and evaluated. + + indices_by_irrep : list or None + Index pairs grouped by irreducible representation when symmetry + is enabled. Otherwise ``None``. + + degens : list or None + Placeholder for degeneracy information. Currently set to ``None`` + when symmetry grouping is not used. + + Notes + ----- + Four index-generation strategies are available: + + ``loop()`` + Generate all upper-triangular coordinate pairs without symmetry. + + ``loop_diagonal()`` + Generate only diagonal coordinate pairs without symmetry. + + ``loop_symmetry()`` + Generate upper-triangular coordinate pairs within each irreducible + representation. + + ``loop_symmetry_diagonal()`` + Generate only diagonal coordinate pairs within each irreducible + representation. + + The :meth:`run` method automatically selects the appropriate strategy + based on the CMA level, derivative level, and symmetry settings. """ def __init__(self, num_deg_free, cma_level, options, proj_irreps=None): @@ -18,6 +92,36 @@ def __init__(self, num_deg_free, cma_level, options, proj_irreps=None): self.proj_irreps = proj_irreps def run(self): + """ + Generate displacement indices for the current CMA calculation. + + Selects and executes the appropriate index-generation algorithm based + on the symmetry settings, CMA level, and derivative level specified + in the options object. + + Selection logic + --------------- + Symmetry enabled + * Level B, energy finite differences: + :meth:`loop_symmetry` + * Level A or gradient finite differences: + :meth:`loop_symmetry_diagonal` + + Symmetry disabled + * Level B, energy finite differences: + :meth:`loop` + * Level A or gradient finite differences: + :meth:`loop_diagonal` + + Returns + ------- + None + + Notes + ----- + Results are stored in the instance attributes ``indices`` and + ``indices_by_irrep``. + """ if self.options.molsym_symmetry: if self.cma_level == "A" or self.options.deriv_level_b: self.loop_symmetry_diagonal() diff --git a/concordantmodes/cma.py b/concordantmodes/cma.py index 3e5e6ea..c12d39c 100644 --- a/concordantmodes/cma.py +++ b/concordantmodes/cma.py @@ -32,23 +32,92 @@ from concordantmodes.zmat import Zmat -class ConcordantModes(object): +class ConcordantModes: """ - This constant is from: - https://physics.nist.gov/cgi-bin/cuu/Value?hr - If this link dies, find the new link on NIST - for the Hartree to Joule conversion and pop it in there. - There is a standard uncertainty of 0.0000000000085 to the MDYNE_HART constant. - BOHR_ANG: Standard uncertainty of 0.00000000080 + Driver class for the Concordant Modes Algorithm (CMA). + + This class orchestrates the complete vibrational analysis workflow, + including: + + - Reading and processing molecular geometries and internal coordinates. + - Generating symmetry-adapted displacement coordinates via molsym. + - Constructing internal-coordinate Hessians from finite-difference + energies or gradients. + - Performing low-level (Level B) and high-level (Level A) force-constant + calculations within the CMA framework. + - Computing G and F matrices and solving the vibrational GF equations. + - Generating Total Energy Distribution (TED) analyses. + - Transforming force constants between internal and Cartesian coordinates. + - Producing harmonic vibrational frequencies, zero-point vibrational + energies (ZPVE), normal modes, and Molden output files. + + The implementation supports conventional harmonic analyses, projected + coordinate spaces, symmetry-adapted displacements, reduced-displacement + CMA calculations, and optional second-order coordinate transformations. + + Parameters + ---------- + options : object + User-defined options controlling coordinate generation, electronic + structure calculations, finite-difference procedures, symmetry + handling, and output settings. + proj : numpy.ndarray, optional + Projection matrix defining the vibrational coordinate space. An empty + array indicates that the projection will be generated automatically. + extra_indices : list, optional + Additional displacement index pairs used for selective off-diagonal + force-constant evaluation in higher-order CMA procedures. + + Attributes + ---------- + MDYNE_HART : float + Conversion factor from Hartrees to mdyne·Å. + BOHR_ANG : float + Conversion factor from Bohr to Angstrom. + proj : numpy.ndarray + Active projection matrix defining a set of linearly independent internal coordinates. + zmat_obj : Zmat + Molecular internal-coordinate representation. + symm_obj : Symmetry + Symmetry analysis and displacement-reduction handler. + s_vec : SVectors + Internal-coordinate B-matrix and projection-space generator. + TED_obj : TED + Total Energy Distribution analysis object. + disp : TransfDisp + Coordinate displacement and transformation manager. + G : numpy.ndarray + Final transformed G matrix. + F_a : numpy.ndarray + High-level (Level A) force-constant matrix. + F_b : numpy.ndarray + Low-level (Level B) force-constant matrix. + + Notes + ----- + The CMA procedure computes a lower-cost reference Hessian (Level B), + obtains vibrational normal coordinates from a GF analysis, and then + evaluates selected force constants at a higher level of theory (Level A) + in the Level B normal coordinate basis. This approach significantly reduces + the computational cost of high-accuracy vibrational frequency calculations + while retaining much of the accuracy of a full Level A Hessian. """ - def __init__(self, options, proj=None, extra_indices=[]): + def __init__(self, options, proj=np.array([]), extra_indices=[]): self.options = options + # These constants are from: + # https://physics.nist.gov/cgi-bin/cuu/Value?hr + # If this link dies, find the new link on NIST + # for the Hartree to Joule conversion and pop it in there. + # There is a standard uncertainty of 0.0000000000085 to the MDYNE_HART constant. + # BOHR_ANG: Standard uncertainty of 0.00000000080 self.MDYNE_HART = 4.3597447222071 self.BOHR_ANG = 0.529177210903 self.proj = proj self.extra_indices = extra_indices + # Change to 'cma' rather than run. We should then be able to introduce a function for higher order + # force constant computations. def run(self, sym_sort=[]): t1 = time.time() @@ -57,6 +126,8 @@ def run(self, sym_sort=[]): # string value to access "state" of CMA program. Are we in level A, B, C, etc... cma_level = "B" + self.sym_sort = sym_sort + # Parse the output to get all pertinent ZMAT info self.zmat_obj = Zmat(self.options) self.zmat_obj.run() @@ -75,275 +146,53 @@ def run(self, sym_sort=[]): self.symm_obj.dummy_obj() self.symm_obj.symtext = None # check if sym_sort object was passed in. If so, intialize sym_sort objects - if len(sym_sort) > 1: - self.symm_obj.create_flat_sym_sort(sym_sort) + if len(self.sym_sort) > 1: + self.symm_obj.create_flat_sym_sort(self.sym_sort) - # Compute the Level B s-vectors - s_vec = SVectors( - self.zmat_obj, self.options # , self.zmat_obj.variable_dictionary_b - ) - s_vec.run( + coord_type = "internal" + if self.options.cart_fc_b: + coord_type = "cartesian" + self.F_b, self.grad_b = self.compute_hessian( + cma_level, + self.options.deriv_level_b, + coord_type, + self.zmat_obj, + self.options, + rootdir, self.zmat_obj.cartesians_b, - True, proj=self.proj, - second_order=self.options.second_order, + off_diag=0, + prog=self.options.program_b, + gen_disps=self.options.gen_disps_b, + calc=self.options.calc_b, ) - if self.options.molsym_symmetry: - self.symm_obj.make_proj(s_vec) - s_vec.proj = copy.deepcopy(self.symm_obj.salc_proj) - - self.TED_obj = TED(s_vec.proj, self.zmat_obj, self.options) - - # Compute G-Matrix - g_mat = GMatrix(self.zmat_obj, s_vec, self.options) + g_mat = GMatrix(self.zmat_obj, self.s_vec, self.options, proj=self.proj) g_mat.run() - G = g_mat.G.copy() - if os.path.exists(rootdir + "/fc_b.grad"): - g_read_obj = GrRead("fc_b.grad") - g_read_obj.run(self.zmat_obj.cartesians_b) - - num_deg_free = s_vec.proj.shape[1] - # Read in FC matrix in cartesians, then convert to internals. - # Or compute an initial hessian in internal coordinates. - self.options.init_bool = False - if os.path.exists(rootdir + "/fc_b.dat"): - f_read_obj = FcRead("fc_b.dat") - self.options.cart_fc_b = True - elif os.path.exists(rootdir + "/FCMFINAL"): - f_read_obj = FcRead("FCMFINAL") - else: - self.options.init_bool = True - - # First generate displacements in internal coordinates - eigs_b = np.eye(len(s_vec.proj.T)) - coord_type = "internal" - if self.options.cart_fc_b: - eigs_b = np.eye(len(self.zmat_obj.cartesians_b.flatten())) - coord_type = "cartesian" - - algo = Algorithm( - num_deg_free, - cma_level, - self.options, - proj_irreps=self.symm_obj.proj_irreps, - ) - algo.run() - if not self.options.deriv_level_b and not self.options.molsym_symmetry: - # Sym_sort doesn't seem to be working - print("symmetric displacements:") - if len(sym_sort) > 1: - algo.indices = self.symm_obj.create_sym_sort_disps( - sym_sort, algo.indices - ) - else: - self.symm_obj.indices_by_irrep = algo.indices_by_irrep - - b_disp = TransfDisp( - s_vec, - self.zmat_obj, - eigs_b, - True, - self.TED_obj, - self.options, - algo.indices, - symm_obj=self.symm_obj, - coord_type=coord_type, - deriv_level=self.options.deriv_level_b, - cma_level=cma_level, - ) - b_disp.run() - - prog_b = self.options.program_b - prog_name_b = prog_b.split("@")[0] - - if self.options.gen_disps_b: - - ref_geom_b = b_disp.disp_cart["ref"] - - dir_obj_b = DirectoryTree( - prog_name_b, - self.zmat_obj, - ref_geom_b, - cma_level, - b_disp.p_disp, - b_disp.m_disp, - self.options, - algo.indices, - self.symm_obj, - "templateB.dat", - "DispsB", - deriv_level=self.options.deriv_level_b, - ) - dir_obj_b.run() - - if not self.options.calc_b: - print( - "The Level B displacements have been generated, now they must be run locally." - ) - raise RuntimeError - else: - if not os.path.exists(rootdir + "/DispsB"): - print( - "You need to have a DispsB directory already present if you want to proceed under current conditions!" - ) - raise RuntimeError - - os.chdir(rootdir + "/DispsB") - - if self.options.calc_b: - sub = Submit(self.options, cma_level, rootdir, prog_name_b, prog_b) - sub.run() - - reap_obj_b = Reap( - self.options, - len(eigs_b), - algo.indices, - self.symm_obj, - cma_level, - deriv_level=self.options.deriv_level_b, - disp=b_disp, - ted=self.TED_obj, - zmat=self.zmat_obj, - ) - reap_obj_b.run() - - p_array_b = reap_obj_b.p_en_array - m_array_b = reap_obj_b.m_en_array - # can this be folded into reap obj? - if not self.options.deriv_level_b: - ref_en_b = reap_obj_b.ref_en - else: - cart_p_array_b = reap_obj_b.p_grad_array - cart_m_array_b = reap_obj_b.m_grad_array - p_array_b = np.zeros(np.eye(len(eigs_b)).shape) - m_array_b = np.zeros(np.eye(len(eigs_b)).shape) - ref_en_b = None - - # Need to convert this array here from cartesians to internals using projected A-tensor - # Surely this can be folded into reap_obj. Then we simply define p_array and m_array so - # so force_constant can read them in and determine what do do with them based off of - # deriv_level. - for i in range(len(algo.indices)): - grad_s_vec = SVectors( - self.zmat_obj, - self.options, - ) - - grad_s_vec.run(b_disp.p_disp[i], False) - A_proj = np.dot(LA.pinv(grad_s_vec.B), self.TED_obj.proj) - p_array_b[i] = np.dot(cart_p_array_b[i].T, A_proj) - - grad_s_vec.run(b_disp.m_disp[i], False) - A_proj = np.dot(LA.pinv(grad_s_vec.B), self.TED_obj.proj) - m_array_b[i] = np.dot(cart_m_array_b[i].T, A_proj) - - fc_b = ForceConstant( - b_disp, - p_array_b, - m_array_b, - ref_en_b, - self.options, - algo.indices, - deriv_level=self.options.deriv_level_b, - coord_type_b=coord_type, - cma_level=cma_level, - ) - fc_b.run() - - f_conv_test = FcConv( - fc_b.FC, - s_vec, - self.zmat_obj, - "cartesian", - False, - self.TED_obj, - self.options, - ) - f_conv_test.run(grad=fc_b.gradient) - - if self.options.second_order and self.options.cart_fc_b: - f_conv_obj = FcConv( - fc_b.FC, - s_vec, - self.zmat_obj, - "internal", - False, - self.TED_obj, - self.options, - ) - f_conv_obj.run(grad=fc_b.gradient) - fc_b.FC = f_conv_obj.F - - if not self.options.init_bool: - f_read_obj.run() - f_conv_obj = FcConv( - f_read_obj.fc_mat, - s_vec, - self.zmat_obj, - "internal", - False, - self.TED_obj, - self.options, - ) - if self.options.second_order: - f_conv_obj.run(grad=g_read_obj.cart_grad) - else: - f_conv_obj.run() - F = f_conv_obj.F - else: - F = fc_b.FC - - # Fold into g_mat - if self.options.coords != "ZMAT": - g_mat.G = np.dot(self.TED_obj.proj.T, np.dot(g_mat.G, self.TED_obj.proj)) - - # These can be folded into prior functions self.options.init_bool = False - F[np.abs(F) < self.options.tol] = 0 - g_mat.G[np.abs(g_mat.G) < self.options.tol] = 0 - del_tol = 1.0e-3 - for row in g_mat.G: - abs_row = np.abs(row) - row[abs_row < np.max(abs_row) * del_tol] = 0 - g_mat.G = g_mat.G.T - for row in g_mat.G: - abs_row = np.abs(row) - row[abs_row < np.max(abs_row) * del_tol] = 0 - g_mat.G = g_mat.G.T - - for row in F: - abs_row = np.abs(row) - row[abs_row < np.max(abs_row) * del_tol] = 0 - F = F.T - for row in F: - abs_row = np.abs(row) - row[abs_row < np.max(abs_row) * del_tol] = 0 - F = F.T - # raise RuntimeError - if len(sym_sort) > 1: - F, g_mat.G = self.symm_obj.GF_sym_sort(F, g_mat, sym_sort) + if len(self.sym_sort) > 1: + F, g_mat.G = self.symm_obj.GF_sym_sort(self.F_b, g_mat, self.sym_sort) # Run the GF matrix method with the internal F-Matrix and computed G-Matrix! print("Level B Frequencies:") b_GF = GFMethod( g_mat.G.copy(), - F.copy(), + self.F_b.copy(), self.zmat_obj, self.TED_obj, self.options, - self.symm_obj.symtext, + symtext=self.symm_obj.symtext, ) b_GF.run() ted_b = b_GF.ted.TED # more sym_sort stuff - if len(sym_sort): + if len(self.sym_sort): self.irreps_b, flat_sym_freqs = self.symm_obj.mode_symmetry_sort( - b_GF.ted.TED, sym_sort, b_GF.freq + ted_b, self.sym_sort, b_GF.freq ) self.ref_b = np.array(flat_sym_freqs) #### this block could probably be moved inside the symmetry.py module? @@ -356,11 +205,12 @@ def run(self, sym_sort=[]): ted_b = ted_b.T #### end of block that could probably be moved inside the symmetry.py module? + self.F_b = np.dot(np.dot(b_GF.L.T, self.F_b), b_GF.L) # Now for the TED check. if self.options.ted_check: self.G = np.dot(np.dot(LA.inv(b_GF.L), g_mat.G), LA.inv(b_GF.L).T) self.G[np.abs(self.G) < self.options.tol] = 0 - self.F = np.dot(np.dot(b_GF.L.T, F), b_GF.L) + self.F = np.dot(np.dot(b_GF.L.T, self.F_b), b_GF.L) self.F[np.abs(self.F) < self.options.tol] = 0 print("TED Frequencies: Degeneracy x Irrep") @@ -403,7 +253,7 @@ def run(self, sym_sort=[]): self.options, algo.indices, deriv_level=self.options.deriv_level_b, - coord_type_b=coord_type, + coord_type=coord_type, cma_level=cma_level, ) FC_c.run() @@ -411,9 +261,8 @@ def run(self, sym_sort=[]): fc_c = FC_c.FC fc_c = np.dot(np.dot(b_GF.L.T, fc_c), b_GF.L) fc_c[np.abs(fc_c) < self.options.tol] = 0 - # print("Level C normal mode Force Constants:") - # print(fc_c) xi = copy.deepcopy(fc_c) * 0.0 + print("xi values:") for i in range(len(xi)): for j in range(i + 1): if i != j: @@ -422,144 +271,51 @@ def run(self, sym_sort=[]): xi[i, j] = buff / np.sqrt( np.abs(fc_c[i, i]) * np.abs(fc_c[j, j]) ) + print(i, j) + print(xi[i, j]) if xi[i, j] > self.options.xi_tol: self.extra_indices.append([j, i]) print("CMA-2 extra off-diag indices:") + print(len(self.extra_indices)) print(self.extra_indices) + # elif self.options.off_diag == '1': + # self.extra_indices = self.od_inds + # Can we generalize compute_hessian to run this too? # Now switch state to cma_level = "A" cma_level = "A" - if self.options.molsym_symmetry: - algo = Algorithm( - num_deg_free, cma_level, self.options, self.symm_obj.proj_irreps - ) - else: - algo = Algorithm(num_deg_free, cma_level, self.options, None) - - algo.run() - if self.options.molsym_symmetry: - self.symm_obj.indices_by_irrep = algo.indices_by_irrep - if len(self.extra_indices): - algo.indices = np.append(algo.indices, self.extra_indices, axis=0) - # Recompute the B-Tensors to match the final geometry, - # then generate the displacements. - - s_vec = SVectors( - self.zmat_obj, self.options # , self.zmat_obj.variable_dictionary_a - ) - s_vec.run( - self.zmat_obj.cartesians_a, - False, - proj=self.TED_obj.proj, - second_order=self.options.second_order, - ) - - transf_disp = TransfDisp( - s_vec, - self.zmat_obj, - b_GF.L, - True, - self.TED_obj, - self.options, - algo.indices, - symm_obj=self.symm_obj, - cma_level=cma_level, - ) - transf_disp.run(fc=F) - p_disp = transf_disp.p_disp - m_disp = transf_disp.m_disp - if self.options.disp_check: - raise RuntimeError - - # The displacements have been generated, now we have to run them! - prog_a = self.options.program_a - prog_name_a = prog_a.split("@")[0] - if self.options.gen_disps_a: - - ref_geom_a = transf_disp.disp_cart["ref"] - dir_obj_a = DirectoryTree( - prog_name_a, - self.zmat_obj, - ref_geom_a, - cma_level, - p_disp, - m_disp, - self.options, - algo.indices, - self.symm_obj, - "templateA.dat", - "DispsA", - ) - dir_obj_a.run() - - if not self.options.calc_a: - print( - "The displacements have been generated, now they must be run locally." - ) - raise RuntimeError - else: - if not os.path.exists(rootdir + "/DispsA"): - print( - "You need to have a DispsA directory already present if you want to proceed under current conditions!" - ) - raise RuntimeError - - os.chdir(rootdir + "/DispsA") - - if self.options.calc_a: - sub = Submit(self.options, cma_level, rootdir, prog_name_a, prog_a) - sub.run() - - # After this point, all of the jobs have finished, and it's time - # to reap the energies as well as checking for sucesses - reap_obj_a = Reap( - self.options, - num_deg_free, - algo.indices, - self.symm_obj, + fc = np.array([]) + if self.options.reduced_disp: + fc = self.F_b + self.F_a, self.grad_a = self.compute_hessian( cma_level, - ) - reap_obj_a.run() - - p_en_array_a = reap_obj_a.p_en_array - m_en_array_a = reap_obj_a.m_en_array - ref_en_a = reap_obj_a.ref_en - - fc_a = ForceConstant( - transf_disp, - p_en_array_a, - m_en_array_a, - ref_en_a, + self.options.deriv_level_a, + "internal", + self.zmat_obj, self.options, - algo.indices, - cma_level=cma_level, + rootdir, + self.zmat_obj.cartesians_a, + proj=self.s_vec.proj, + off_diag=self.options.off_diag, + prog=self.options.program_a, + gen_disps=self.options.gen_disps_a, + calc=self.options.calc_a, + eigs=b_GF.L, + fc=fc, ) - fc_a.run() - - np.set_printoptions(precision=4, linewidth=240) - # print("Computed Force Constants:") - # print(fc_a.FC) - - # print("Computed Normal Mode Gradient:") - # print(fc_a.gradient) - - self.F = fc_a.FC # Recompute the G-matrix with the new geometry, and then transform # the G-matrix using the lower level of theory eigenvalue matrix. # This will not fully diagonalize the G-matrix if a different # geometry is used between the two. - - g_mat = GMatrix(self.zmat_obj, s_vec, self.options) + g_mat = GMatrix(self.zmat_obj, self.s_vec, self.options, proj=self.proj) g_mat.run() - - if self.options.coords != "ZMAT": - g_mat.G = np.dot(self.TED_obj.proj.T, np.dot(g_mat.G, self.TED_obj.proj)) + np.set_printoptions(precision=7, linewidth=240) self.G = g_mat.G - self.G = np.dot(np.dot(transf_disp.eig_inv, self.G), transf_disp.eig_inv.T) + self.G = np.dot(np.dot(self.disp.eig_inv, self.G), self.disp.eig_inv.T) self.G[np.abs(self.G) < self.options.tol] = 0 if self.options.benchmark_full: @@ -571,7 +327,7 @@ def run(self, sym_sort=[]): print("Final Harmonic Frequencies:") a_GF = GFMethod( self.G, - self.F, + self.F_a, self.zmat_obj, self.TED_obj, self.options, @@ -603,16 +359,16 @@ def run(self, sym_sort=[]): # coordinates and writes out an "output.default.hess" file, which # is of the same format as FCMFINAL of CFOUR. - self.F = np.dot(np.dot(transf_disp.eig_inv.T, self.F), transf_disp.eig_inv) - self.gradient = np.dot(fc_a.gradient, transf_disp.eig_inv) + self.F_a = np.dot(np.dot(self.disp.eig_inv.T, self.F_a), self.disp.eig_inv) + self.gradient = np.dot(self.grad_a, self.disp.eig_inv) cart_conv = FcConv( - self.F, - s_vec, + self.F_a, + self.s_vec, self.zmat_obj, "cartesian", True, - self.TED_obj, + self.proj, self.options, ) if self.options.second_order: @@ -623,19 +379,245 @@ def run(self, sym_sort=[]): self.F_cart = cart_conv.F - # if mol2.size != 0: - # if self.options.rmsd: - # rmsd_geom = RMSD() - # rmsd_geom.run(mol1,mol2) - print("Frequency Shift (cm^-1): ") print(a_GF.freq - b_GF.freq) for i in a_GF.freq - b_GF.freq: print(i) + print("RMSD Freq Shift (cm^-1): ") + print(np.sqrt(np.mean((a_GF.freq - b_GF.freq) ** 2))) + print("MAX Freq Shift (cm^-1): ") + print(np.max(np.abs(a_GF.freq - b_GF.freq))) # Write a molden file - molden = MoldenWriter(self.zmat_obj, transf_disp, a_GF.freq) + molden = MoldenWriter(self.zmat_obj, self.disp, a_GF.freq) molden.run() t2 = time.time() print("This program took " + str(t2 - t1) + " seconds to run.") + + # It's time to make this a real driver. Begin by creating a procedure for computing an arbitrary hessian + # whether it's a full hessian, variable derivative levels, variable CMA_levels, variable coordinate types. + def compute_hessian( + self, + cma_level, + deriv_level, + coord_type, + zmat, + options, + rootdir, + ref_carts, + proj=np.array([]), + off_diag=None, + prog=None, + gen_disps=True, + calc=True, + eigs=None, + fc=np.array([]), + ): + self.s_vec = SVectors(zmat, options) + b_proj = False + if cma_level == "B" and not options.man_proj: + b_proj = True + self.s_vec.run( + ref_carts, + b_proj, + proj=proj, + second_order=options.second_order, + ) + + if self.options.molsym_symmetry: + self.symm_obj.make_proj(self.s_vec) + self.s_vec.proj = copy.deepcopy(self.symm_obj.salc_proj) + + self.proj = self.s_vec.proj + + self.TED_obj = TED(self.proj, zmat, options) + + if os.path.exists(rootdir + "/fc_" + cma_level.lower() + ".grad"): + g_read_obj = GrRead("fc_" + cma_level.lower() + ".grad") + # Need to pass in general carts here. + g_read_obj.run(ref_carts) + + num_deg_free = self.proj.shape[1] + if deriv_level and self.options.cart_fc_b: + num_deg_free = len(ref_carts.flatten()) + options.init_bool = False + cart_fc = False + if os.path.exists(rootdir + "/fc_" + cma_level.lower() + ".dat"): + f_read_obj = FcRead("fc_" + cma_level.lower() + ".dat") + cart_fc = True + + f_read_obj.run() + f_conv_obj = FcConv( + f_read_obj.fc_mat, + self.s_vec, + zmat, + "internal", + False, + self.TED_obj, + options, + ) + if self.options.second_order: + f_conv_obj.run(grad=g_read_obj.cart_grad) + else: + f_conv_obj.run() + F = f_conv_obj.F + else: + options.init_bool = True + + # First generate displacements in internal coordinates + if cma_level == "B": + eigs = np.eye(len(self.proj.T)) + if self.options.cart_fc_b: + eigs = np.eye(len(ref_carts.flatten())) + cart_fc = True + + algo = Algorithm( + num_deg_free, + cma_level, + options, + proj_irreps=self.symm_obj.proj_irreps, + ) + algo.run() + print("Pre sym indices:") + print(algo.indices) + # raise RuntimeError + if ( + not options.deriv_level_b + and not options.molsym_symmetry + and cma_level == "B" + ): + # Sym_sort doesn't seem to be working + print("symmetric displacements:") + if len(self.sym_sort) > 1: + algo.indices = self.symm_obj.create_sym_sort_disps( + self.sym_sort, algo.indices + ) + else: + self.symm_obj.indices_by_irrep = algo.indices_by_irrep + print("Post sym indices:") + print(algo.indices) + if cma_level == "A" and len(self.extra_indices): + algo.indices += self.extra_indices + self.disp = TransfDisp( + self.s_vec, + zmat, + eigs, + True, + self.proj, + options, + algo.indices, + symm_obj=self.symm_obj, + coord_type=coord_type, + deriv_level=deriv_level, + cma_level=cma_level, + ) + self.disp.run(fc=fc) + + prog_name = prog.split("@")[0] + + if gen_disps: + ref_geom = self.disp.disp_cart["ref"] + + dir_obj = DirectoryTree( + prog_name, + zmat, + ref_geom, + cma_level, + self.disp.p_disp, + self.disp.m_disp, + options, + algo.indices, + self.symm_obj, + "template" + cma_level.upper() + ".dat", + "Disps" + cma_level.upper(), + deriv_level=deriv_level, + ) + dir_obj.run() + + if not calc: + print( + # "The Level B displacements have been generated, now they must be run locally." + "The Level " + + cma_level.upper() + + " displacements have been generated, now they must be run locally." + ) + raise RuntimeError + else: + if not os.path.exists(rootdir + "/Disps" + cma_level.upper()): + print( + "You need to have a Disps" + + cma_level.upper() + + " directory already present if you want to proceed under current conditions!" + ) + raise RuntimeError + + if calc: + sub = Submit(options, cma_level, rootdir, prog_name, prog) + sub.run() + + os.chdir(rootdir + "/Disps" + cma_level.upper()) + + reap_obj = Reap( + options, + len(eigs), + algo.indices, + self.symm_obj, + cma_level, + deriv_level=deriv_level, + disp=self.disp, + proj=self.proj, + zmat=zmat, + ) + reap_obj.run() + os.chdir(rootdir) + # can this be folded into reap obj? + if not deriv_level: + ref_en = reap_obj.ref_en + p_array = reap_obj.p_en_array + m_array = reap_obj.m_en_array + else: + p_array = reap_obj.p_grad_array + m_array = reap_obj.m_grad_array + ref_en = None + + fc = ForceConstant( + self.disp, + p_array, + m_array, + ref_en, + options, + algo.indices, + deriv_level=deriv_level, + coord_type=coord_type, + cma_level=cma_level, + gradient=reap_obj.ref_grad, + ) + fc.run() + + if options.second_order and cart_fc: + f_conv_obj = FcConv( + fc.FC, + self.s_vec, + zmat, + "internal", + False, + self.proj, + options, + ) + f_conv_obj.run(grad=fc.gradient) + fc.FC = f_conv_obj.F + + F = fc.FC + + F[np.abs(F) < self.options.tol] = 0 + del_tol = 1.0e-3 + for row in F: + abs_row = np.abs(row) + row[abs_row < np.max(abs_row) * del_tol] = 0 + F = F.T + for row in F: + abs_row = np.abs(row) + row[abs_row < np.max(abs_row) * del_tol] = 0 + F = F.T + return F, fc.gradient diff --git a/concordantmodes/directory_tree.py b/concordantmodes/directory_tree.py index 5c4f32c..6927848 100644 --- a/concordantmodes/directory_tree.py +++ b/concordantmodes/directory_tree.py @@ -1,12 +1,137 @@ -import fileinput -import os -import re +from pathlib import Path import shutil import copy -import numpy as np -class DirectoryTree(object): +class DirectoryTree: + """ + Generate and organize displacement calculation directories for the + Concordant Modes Algorithm (CMA). + + This class creates the directory structure and electronic structure + program input files required for finite-difference force-constant + calculations. For each reference geometry and displaced geometry, a + numbered subdirectory is created containing an appropriately formatted + input file and any required supporting files. + + Supported quantum chemistry programs include Molpro, Psi4, CFOUR, + and ORCA. Geometry coordinates are inserted into a user-provided + template file at a configurable location. + + Depending on the derivative level, the class generates inputs for: + + - Energy finite differences (Hessians from energies) + - Gradient finite differences (Hessians from gradients) + + The generated directory structure is compatible with the CMA workflow + and subsequent job-submission and result-harvesting stages. + + Parameters + ---------- + prog_name : str + Electronic structure program name. Supported values are + ``"molpro"``, ``"psi4"``, ``"cfour"``, and ``"orca"``. + + zmat : Zmat + Molecular coordinate object containing atom labels and geometry + information. + + ref_geom : ndarray of shape (n_atoms, 3) + Reference Cartesian geometry. + + cma_level : {"A", "B", "C"} + CMA stage for which displacement calculations are being generated. + + - ``"A"`` : High-level force constants. + - ``"B"`` : Low-level reference force constants. + - ``"C"`` : Optional off-diagonal refinement calculations. + + p_disp : ndarray or dict + Positive displacement geometries. + + m_disp : ndarray or dict + Negative displacement geometries. + + options : object + User options object containing template insertion locations and + symmetry settings. + + indices : list[list[int]] + Displacement coordinate indices generated by the CMA algorithm. + + symm_obj : Symmetry + Symmetry analysis object used to determine symmetry-reduced + displacement sets. + + template : str or pathlib.Path + Input template file containing all program keywords and settings. + Cartesian coordinates are inserted at the specified location. + + dir_name : str or pathlib.Path + Name of the root directory in which displacement calculations will + be generated. + + deriv_level : int, optional + Finite-difference derivative level. + + - ``0`` : Energy finite differences. + - ``1`` : Gradient finite differences. + + Attributes + ---------- + insertion_index : int + Line number in the template where Cartesian coordinates are + inserted. + + dir_name : pathlib.Path + Root displacement directory. + + template : pathlib.Path + Template input file. + + init : bool + Indicates whether an ``initden.dat`` file is available. + + genbas : bool + Indicates whether a ``GENBAS`` file is available. + + ecp : bool + Indicates whether an ``ECPDATA`` file is available. + + sub : bool + Indicates whether a ``sub_script.sh`` file is available. + + Notes + ----- + Directory numbering follows the CMA convention: + + - Directory ``1`` contains the reference geometry. + - Subsequent directories contain positive and negative displacements. + - For energy finite differences, displacement pairs correspond to + force-constant index pairs ``(i,j)``. + - For gradient finite differences, displacement pairs correspond to + individual coordinates. + + If an existing displacement directory is present, it is copied to an + ``old`` backup before a new directory tree is created. + + Optional files detected in the working directory are automatically + copied into every displacement directory: + + - ``initden.dat`` + - ``GENBAS`` + - ``ECPDATA`` + - ``sub_script.sh`` + """ + + PROG_LIST = {"molpro", "psi4", "cfour", "orca"} + + INSERTION_MAP = { + "A": "cart_insert_a", + "B": "cart_insert_b", + "C": "cart_insert_c", + } + def __init__( self, prog_name, @@ -31,206 +156,136 @@ def __init__( self.options = options self.indices = indices self.symm_obj = symm_obj - self.template = template - self.dir_name = dir_name + self.template = Path(template) + self.dir_name = Path(dir_name) self.deriv_level = deriv_level - if cma_level == "A": - self.insertion_index = self.options.cart_insert_a - elif cma_level == "B": - self.insertion_index = self.options.cart_insert_b - elif cma_level == "C": - self.insertion_index = self.options.cart_insert_c - else: - print("Please specify A, B, or C for your cma_level") - raise RuntimeError + + try: + self.insertion_index = getattr(self.options, self.INSERTION_MAP[cma_level]) + except KeyError: + raise RuntimeError("cma_level must be A, B, or C") def run(self): + """ + Create the complete displacement directory tree. + + Reads the template input file, validates the selected electronic + structure program, backs up any existing displacement directory, + determines available support files, and generates all reference and + displaced-geometry input files required for the requested derivative + level. - root = os.getcwd() + Returns + ------- + None - prog_name = self.prog_name + Raises + ------ + RuntimeError + If an unsupported quantum chemistry program is specified. + RuntimeError + If an unsupported derivative level is requested. + """ + root = Path.cwd() + + if self.prog_name not in self.PROG_LIST: + raise RuntimeError(f"Unsupported program: {self.prog_name}") + + data = self.template.read_text().splitlines(keepends=True) n_atoms = len(self.zmat.atom_list) - prog_list = ["molpro", "psi4", "cfour", "orca"] + # Detect optional files + self.init = (root / "initden.dat").exists() + self.genbas = (root / "GENBAS").exists() + self.ecp = (root / "ECPDATA").exists() + self.sub = (root / "sub_script.sh").exists() - if prog_name in prog_list: - with open(self.template, "r") as file: - data = file.readlines() - else: - print("Specified program not supported: " + prog_name) - raise RuntimeError - - self.init = False - self.genbas = False - self.ecp = False - self.sub = False - print("Giraffe") - print(root) - if os.path.exists(root + "/initden.dat"): - self.init = True - if os.path.exists(root + "/GENBAS"): - self.genbas = True - if os.path.exists(root + "/ECPDATA"): - self.ecp = True - if os.path.exists(root + "/sub_script.sh"): - print("We made it here.") - self.sub = True - # raise RuntimeError - - data_buff = data.copy() - if os.path.exists(os.getcwd() + "/old" + self.dir_name): - shutil.rmtree("old" + self.dir_name, ignore_errors=True) - if os.path.exists(os.getcwd() + "/" + self.dir_name): - shutil.copytree( - os.getcwd() + "/" + self.dir_name, os.getcwd() + "/old" + self.dir_name - ) - shutil.rmtree(os.getcwd() + "/" + self.dir_name) - inp = "" - if self.prog_name == "cfour": - inp = "ZMAT" - else: - inp = "input.dat" - if os.path.exists(os.getcwd() + "/" + self.dir_name): - shutil.move(self.dir_name, "old" + self.dir_name) - os.mkdir(self.dir_name) - os.chdir("./" + self.dir_name) - if not self.deriv_level: - self.make_input( - copy.deepcopy(data), - self.ref_geom, - n_atoms, - self.zmat.atom_list, - self.insertion_index, - inp, - "1", - ) - - # Not including the reference input, this function generates the directories for the displacement - # jobs and copies in the input file data. Following this, these jobs are ready to be submitted to the queue. - - p_disp = self.p_disp - m_disp = self.m_disp - indices = self.indices - direc = 2 - - if self.symm_obj.symtext is not None and self.options.exploit_pm_symm: - if self.options.only_TSIR: - indices = self.symm_obj.indices_by_irrep[0] - else: - indices = ( - deepcopy.copy(self.symm_obj.indices_by_irrep) - .flatten() - .reshape((-1, 2)) - ) - - Sum = 1 - h = 0 - for index in indices: - i, j = index[0], index[1] - self.make_input( - copy.deepcopy(data), - p_disp[i, j], - n_atoms, - self.zmat.atom_list, - self.insertion_index, - inp, - direc, - ) - if h != 0: - direc += 1 - else: - self.make_input( - copy.deepcopy(data), - m_disp[i, j], - n_atoms, - self.zmat.atom_list, - self.insertion_index, - inp, - direc + 1, - ) - - direc += 2 - if ( - self.symm_obj.symtext is not None - and self.options.exploit_pm_symm - and not self.options.only_TSIR - and Sum > len(self.symm_obj.indices_by_irrep[0]) - ): - h += 1 - Sum += 1 + # Backup existing directory + old_dir = root / f"old{self.dir_name}" + if old_dir.exists(): + shutil.rmtree(old_dir) + if self.dir_name.exists(): + shutil.copytree(self.dir_name, old_dir) + shutil.rmtree(self.dir_name) + + self.dir_name.mkdir() + Path.chdir = lambda p: None # avoid lint complaints + + inp = "ZMAT" if self.prog_name == "cfour" else "input.dat" + + # Generate inputs + if self.deriv_level == 0: + self._run_energy(data, n_atoms, inp) elif self.deriv_level == 1: - direc = 1 - for index in self.indices: - self.make_input( - copy.deepcopy(data), - self.p_disp[index[0]], - n_atoms, - self.zmat.atom_list, - self.insertion_index, - inp, - direc, - ) - - self.make_input( - copy.deepcopy(data), - self.m_disp[index[0]], - n_atoms, - self.zmat.atom_list, - self.insertion_index, - inp, - direc + 1, - ) - - direc += 2 - else: - print( - "Only energy and gradient derivatives are supported. Check your deriv_level keyword." - ) - raise RuntimeError - - def make_input(self, data, dispp, n_at, at, index, inp, direc, copy_files=False): - os.mkdir(str(direc)) - os.chdir("./" + str(direc)) - data_buff = copy.deepcopy(data) - space = " " - if self.prog_name == "cfour": - space = "" - if index == -1: - print( - "The user needs to specify a value for the \ - cart_insert_a, cart_insert_b, or cart_insert_c keyword." - ) - raise RuntimeError + self._run_gradient(data, n_atoms, inp) else: - for i in range(n_at): - data.insert( - index + i, - space - + at[i] - + "{:16.10f}".format(dispp[i][0]) - + "{:16.10f}".format(dispp[i][1]) - + "{:16.10f}".format(dispp[i][2]) - + "\n", - ) - - with open(inp, "w") as file: - file.writelines(data) - data_final = copy.deepcopy(data) - data = copy.deepcopy(data_buff) - if self.options.cluster.lower() == "custom": - copy_files = True - if copy_files: - if self.init: - shutil.copy("../../initden.dat", ".") - if self.genbas: - shutil.copy("../../GENBAS", ".") - if self.ecp: - shutil.copy("../../ECPDATA", ".") - if self.sub: - shutil.copy("../../sub_script.sh", ".") - print("we also made it here") - os.chdir("..") - - return data_final + raise RuntimeError("Only energy and gradient derivatives supported") + + def _run_energy(self, data, n_atoms, inp): + self.make_input(data, self.ref_geom, n_atoms, inp, "1") + + indices = self._resolve_indices() + + direc = 2 + for i, j in indices: + self.make_input(data, self.p_disp[i, j], n_atoms, inp, str(direc)) + self.make_input(data, self.m_disp[i, j], n_atoms, inp, str(direc + 1)) + direc += 2 + + def _run_gradient(self, data, n_atoms, inp): + + self.make_input(data, self.ref_geom, n_atoms, inp, "1") + + direc = 2 + + for idx in self.indices: + self.make_input(data, self.p_disp[idx[0]], n_atoms, inp, str(direc)) + self.make_input(data, self.m_disp[idx[0]], n_atoms, inp, str(direc + 1)) + direc += 2 + + def _resolve_indices(self): + if self.symm_obj.symtext is not None and self.options.exploit_pm_symm: + if self.options.only_TSIR: + return self.symm_obj.indices_by_irrep[0] + return self.symm_obj.indices_by_irrep.reshape((-1, 2)) + + return self.indices + + def make_input(self, data, geom, n_atoms, inp, direc): + path = self.dir_name / str(direc) + path.mkdir(parents=True, exist_ok=True) + + space = "" if self.prog_name == "cfour" else " " + + new_data = data.copy() + + if self.insertion_index == -1: + raise RuntimeError("Invalid insertion index") + + for i in range(n_atoms): + atom = self.zmat.atom_list[i] + x, y, z = geom[i] + line = f"{space}{atom}{x:16.10f}{y:16.10f}{z:16.10f}\n" + new_data.insert(self.insertion_index + i, line) + + (path / inp).write_text("".join(new_data)) + + self.new_data = new_data + + self._copy_support_files(path) + + def _copy_support_files(self, path): + root = Path.cwd().parent + + files = { + "initden.dat": self.init, + "GENBAS": self.genbas, + "ECPDATA": self.ecp, + "sub_script.sh": self.sub, + } + + for fname, exists in files.items(): + if exists: + shutil.copy(root / fname, path) diff --git a/concordantmodes/f_convert.py b/concordantmodes/f_convert.py index 4cc2a58..c8a8bf8 100644 --- a/concordantmodes/f_convert.py +++ b/concordantmodes/f_convert.py @@ -3,25 +3,135 @@ from numpy import linalg as LA -class FcConv(object): +class FcConv: """ - This class may be used to convert a cartesian F-matrix to an internal - F-Matrix, and vice versa. - This constant is from: - https://physics.nist.gov/cgi-bin/cuu/Value?hr - If this link dies, find the new link on NIST - for the Hartree to Joule conversion and pop - it in there. - MDYNE_HART: Standard uncertainty of 0.0000000000085 - BOHR_ANG: Standard uncertainty of 0.00000000080 + Transform force-constant matrices between Cartesian and internal + coordinate representations. + + This class performs the coordinate transformation of harmonic force + constants using the Wilson B-matrix formalism. Transformations may be + performed in either direction: + + - Cartesian → Internal coordinates + - Internal → Cartesian coordinates + + The class also supports projection into a reduced vibrational coordinate + space and optional second-order corrections required for non-stationary + geometries when Cartesian Hessians are transformed into internal + coordinates. + + For Cartesian-to-internal transformations, the generalized inverse of + the Wilson B-matrix is constructed using + + Aᵀ = G⁻¹B + + where + + G = BBᵀ + + is the Wilson G-matrix. + + Optional gradient-dependent corrections based on second derivatives of + the internal coordinates (B² tensors) can be included when transforming + Hessians obtained at non-stationary points on the potential energy + surface. + + Parameters + ---------- + fc_mat : ndarray + Input force-constant matrix. The interpretation depends on the + selected coordinate system: + + - Cartesian Hessian when ``coord="internal"`` + - Internal-coordinate force constant matrix when + ``coord="cartesian"`` + + s_vec : SVectors + Internal-coordinate object containing Wilson B-matrices and, + optionally, second-order B tensors. + + zmat : Zmat + Molecular coordinate representation. + + coord : {"internal", "cartesian"} + Desired output coordinate system. + + - ``"internal"`` transforms Cartesian force constants into + internal coordinates. + - ``"cartesian"`` transforms internal force constants into + Cartesian coordinates. + + print_f : bool + If True, write transformed force constants to disk. + + proj : ndarray + Projection matrix defining a reduced vibrational coordinate basis. + An empty array indicates that the full internal-coordinate basis + should be used. + + options : object + User options controlling units, second-order transformations, + and output generation. + + Attributes + ---------- + F : ndarray + Transformed force-constant matrix. + + A_T : ndarray + Transpose of the generalized inverse transformation matrix used + for Cartesian-to-internal transformations. + + grad : ndarray + Transformed gradient vector when second-order corrections are + applied. + + v_q : ndarray + Internal-coordinate gradient vector used in second-order + corrections. + + MDYNE_HART : float + Conversion factor between Hartree and mdyne·Å. + + BOHR_ANG : float + Conversion factor between Bohr and Angstrom. + + Notes + ----- + The Cartesian-to-internal transformation is performed as + + F_int = Aᵀ F_cart A + + while the reverse transformation is + + F_cart = Bᵀ F_int B + + where B is the Wilson B-matrix and A is its generalized inverse. + + When second-order transformations are enabled, the force constants are + corrected using second derivatives of the internal coordinates stored + in the B² tensor. This correction is required whenever the reference + geometry is not a stationary point and gradient contributions do not + vanish. + + The transformed force constants may optionally be written in a format + compatible with CFOUR-style force-constant files. """ - def __init__(self, fc_mat, s_vec, zmat, coord, print_f, ted, options): + # The constants are from: + # https://physics.nist.gov/cgi-bin/cuu/Value?hr + # If this link dies, find the new link on NIST + # for the Hartree to Joule conversion and pop + # it in there. + # MDYNE_HART: Standard uncertainty of 0.0000000000085 + # BOHR_ANG: Standard uncertainty of 0.00000000080 + + def __init__(self, fc_mat, s_vec, zmat, coord, print_f, proj, options): self.coord = coord self.F = fc_mat self.print_f = print_f self.s_vec = s_vec - self.ted = ted + self.proj = proj self.zmat = zmat self.options = options @@ -29,13 +139,38 @@ def __init__(self, fc_mat, s_vec, zmat, coord, print_f, ted, options): self.BOHR_ANG = 0.529177210903 def run(self, grad=np.array([])): + """ + Perform the force-constant coordinate transformation. + + Depending on the value of ``self.coord``, transforms the input + force-constant matrix between Cartesian and internal coordinate + representations. If gradient information is supplied and second-order + transformations are enabled, gradient-dependent corrections are + included. + + Parameters + ---------- + grad : ndarray, optional + Gradient vector at the reference geometry. Required for + second-order transformations involving non-stationary geometries. + + Returns + ------- + None + + Notes + ----- + Results are stored in the instance attribute ``F``. If second-order + corrections are applied, the transformed gradient is stored in + ``grad`` and the correction term is computed from the second-order B tensor. + """ # First construct the transpose of the A matrix. if self.coord.lower() == "internal": # The cartesian force constants must be in units of Hartree/bohr^2. - if self.ted.proj is None: + if not len(self.proj): B = self.s_vec.B else: - B = np.dot(self.ted.proj.T, self.s_vec.B) + B = np.dot(self.proj.T, self.s_vec.B) G = np.dot(B, B.T) self.A_T = np.dot(LA.inv(G), B) if self.options.units == "MdyneAng": @@ -50,11 +185,12 @@ def run(self, grad=np.array([])): # The basis will eventually need to be projected anyways. np.set_printoptions(precision=6, linewidth=240) self.v_q = np.dot(self.A_T, grad) - if self.ted.proj is None: + if not len(self.proj): B2 = self.s_vec.B2 else: - B2 = np.einsum("rp,pij->rij", self.ted.proj.T, self.s_vec.B2) - C2 = np.einsum("rij,pi,qj->rpq", B2, self.A_T, self.A_T) + B2 = np.einsum("rp,pij->rij", self.proj.T, self.s_vec.B2) + C2 = np.einsum("rij,pi->rpj", B2, self.A_T) + C2 = np.einsum("rpj,qj->rpq", C2, self.A_T) V2 = np.einsum("q,qpr->pr", self.v_q, C2) grad = np.dot(grad, self.A_T.T) @@ -65,20 +201,20 @@ def run(self, grad=np.array([])): self.print_const(fc_name="fc_int.dat", grad=grad) elif self.coord.lower() == "cartesian": - if self.ted.proj is None: + if not len(self.proj): B = self.s_vec.B else: - B = np.dot(self.ted.proj.T, self.s_vec.B) + B = np.dot(self.proj.T, self.s_vec.B) self.F = np.einsum("pi,rj,pr->ij", B, B, self.F) V2 = self.F.copy() * 0 if len(grad) and self.options.second_order: - if self.ted.proj is None: + if not len(self.proj): B2 = self.s_vec.B2 else: - B2 = np.einsum("rp,pij->rij", self.ted.proj.T, self.s_vec.B2) + B2 = np.einsum("rp,pij->rij", self.proj.T, self.s_vec.B2) V2 = np.einsum("rij,r->ij", B2, grad) grad = np.dot(grad, B) @@ -86,17 +222,29 @@ def run(self, grad=np.array([])): self.F += V2 - # print("The Bs:") - # print(self.s_vec.B.shape) - # print(self.s_vec.B) - # print(B.shape) - # print(B) - # raise RuntimeError - if self.print_f: self.print_const(grad=grad) def print_const(self, fc_name="fc_a.dat", grad=np.array([])): + """ + Write transformed force constants and gradients to disk. + + The force-constant matrix is flattened and written in a format + compatible with CFOUR-style force-constant files. If a gradient is + provided, it is written to a separate gradient file. + + Parameters + ---------- + fc_name : str, optional + Output force-constant filename. + + grad : ndarray, optional + Gradient vector to be written alongside the force constants. + + Returns + ------- + None + """ self.N = len(self.F) fc_output = "" fc_output += "{:5d}{:5d}\n".format(len(self.zmat.atom_list), self.N) @@ -119,7 +267,6 @@ def print_const(self, fc_name="fc_a.dat", grad=np.array([])): if len(grad): g_print = grad - # g_print = g_print.flatten() gr_output = "" for i in range(len(g_print) // 3): gr_output += "{:20.10f}".format(g_print[3 * i]) diff --git a/concordantmodes/f_read.py b/concordantmodes/f_read.py index b321f34..c60cddf 100644 --- a/concordantmodes/f_read.py +++ b/concordantmodes/f_read.py @@ -4,11 +4,51 @@ from numpy import linalg as LA -class FcRead(object): +class FcRead: """ - This class may be used to read in force constant information from the - standard 3-column formatted file. Idk what that format is called, but - it's pretty standard. + Read a force-constant matrix from a formatted text file. + + This class parses force constants stored in the standard flattened + three-column format commonly used by quantum chemistry programs and + vibrational analysis utilities. All floating-point values are extracted + from the file, assembled into a one-dimensional array, and reshaped into + a square force-constant matrix. + + The matrix dimension is determined automatically from the total number + of values read: + + N = sqrt(number of force constants) + + Parameters + ---------- + file_name : str + Name of the force-constant file to be read. + + Attributes + ---------- + file_name : str + Input filename. + + fc_regex : re.Pattern + Regular expression used to extract floating-point values from the + file. + + fc_mat : ndarray + Square force-constant matrix reconstructed from the file contents. + + Notes + ----- + The reader assumes that the file contains only the numerical elements + of a square force-constant matrix written in row-major order. The total + number of values must therefore be a perfect square. + + Examples + -------- + >>> reader = FcRead("fc.dat") + >>> reader.run() + >>> F = reader.fc_mat + >>> F.shape + (54, 54) """ def __init__(self, file_name): @@ -16,6 +56,26 @@ def __init__(self, file_name): self.fc_regex = re.compile(r"(-?\d+\.\d+)") def run(self): + """ + Read and reconstruct the force-constant matrix. + + Opens the specified force-constant file, extracts all floating-point + values, and reshapes them into a square NumPy array. + + Returns + ------- + None + + Notes + ----- + The resulting matrix is stored in the attribute ``fc_mat``. + + Raises + ------ + ValueError + If the number of extracted values cannot be reshaped into a + square matrix. + """ with open(self.file_name, "r") as file: fc = file.read() diff --git a/concordantmodes/force_constant.py b/concordantmodes/force_constant.py index 3e469c6..a79e8a1 100644 --- a/concordantmodes/force_constant.py +++ b/concordantmodes/force_constant.py @@ -3,21 +3,124 @@ from numpy import linalg as LA -class ForceConstant(object): - # This script will calculate the force constants of the CMA normal - # modes using numerical differentiation. +class ForceConstant: + """ + Compute force constants by finite differentiation of energies or gradients. - # Symmetric First Derivative - # [f(x+h) - f(x-h)] / 2*h = - # [g_(i_plus) - g_(i_minus)] / 2*disp_size + This class constructs harmonic force-constant matrices in the coordinate + basis defined by the displacement generator. Force constants are obtained + using central finite-difference formulas applied to electronic energies + or gradients evaluated at displaced geometries. - # Diagonal Second Derivative - # [f(x+h) - 2f(x) + f(x-h)] / h^2 = - # [E_(i_plus) - 2*E_(ref) + E_(i_minus)] / disp_size^2 + Two differentiation schemes are supported: - # Off-diagonal Second Derivative - # [f(x+h,y+k) - f(x+h,y) - f(x,y+k) + 2f(x,y) - f(x-h,y) - f(x,y-k) + f(x-h,y-k) / 2*h^2 - # [E_(ij_plus) - E_(i_plus,j) - E_(i,j_plus) + E_(ref) - E_(i_minus,j) - E_(i,j_minus) + E_(i_minus,j_minus) / 2*disp_size^2 + * Energy finite differences (``deriv_level=0``) + - Computes first derivatives (gradients) and second derivatives + (force constants) from displaced energies. + - Includes diagonal and off-diagonal Hessian elements. + + * Gradient finite differences (``deriv_level=1``) + - Computes force constants directly from displaced gradients. + - Produces a symmetric Hessian matrix. + + The resulting force-constant matrix is expressed in the displacement + coordinate basis used by the Concordant Modes Algorithm (CMA), which may + correspond to projected internal coordinates, normal coordinates, or + Cartesian coordinates depending on the calculation setup. + + Parameters + ---------- + disp : TransfDisp + Displacement object containing displacement magnitudes and + coordinate transformation information. + + p_array : ndarray + Energies or gradients evaluated at positive displacements. + + m_array : ndarray + Energies or gradients evaluated at negative displacements. + + ref_en : float or None + Reference energy at the equilibrium geometry. Required for + energy-based finite differences and ignored for gradient-based + calculations. + + options : object + User-defined options controlling the CMA calculation. + + indices : list[list[int]] + Coordinate index pairs defining which force constants are to be + evaluated. + + deriv_level : int, optional + Differentiation level used to construct the Hessian. + + - ``0`` : Hessian from energies. + - ``1`` : Hessian from gradients. + + coord_type : {"internal", "cartesian"}, optional + Coordinate system in which the displacements were generated. + + cma_level : {"A", "B", "C"}, optional + CMA stage associated with the force-constant calculation. + + gradient : ndarray, optional + Reference gradient vector. Primarily used for higher-order + coordinate transformations and non-stationary reference points. + + Attributes + ---------- + FC : ndarray + Computed force-constant matrix. + + gradient : ndarray + Gradient vector computed from central finite differences when + ``deriv_level=0``. + + Notes + ----- + The following central finite-difference formulas are employed. + + First derivative: + + .. math:: + + \\frac{df}{dq_i} + = + \\frac{f(q_i+h)-f(q_i-h)} + {2h} + + Diagonal force constant: + + .. math:: + + \\frac{\\partial^2 E}{\\partial q_i^2} + = + \\frac{E(q_i+h)-2E_0+E(q_i-h)} + {h^2} + + Off-diagonal force constant: + + .. math:: + + \\frac{\\partial^2 E} + {\\partial q_i \\partial q_j} + = + \\frac{ + E(+i,+j)-E(+i)-E(+j) + +2E_0 + -E(-i)-E(-j) + +E(-i,-j)} + {2h_i h_j} + + Symmetry of the Hessian is enforced after construction: + + .. math:: + + F_{ij} = F_{ji} + + for all coordinate pairs. + """ def __init__( self, @@ -28,8 +131,9 @@ def __init__( options, indices, deriv_level=0, - coord_type_b="internal", + coord_type="internal", cma_level="B", + gradient=[], ): self.options = options self.disp = disp @@ -38,13 +142,34 @@ def __init__( self.ref_en = ref_en self.indices = indices self.deriv_level = deriv_level - self.coord_type_b = coord_type_b + self.coord_type = coord_type self.cma_level = cma_level + self.gradient = gradient def run(self): + """ + Construct the force-constant matrix. + + Depending on the selected derivative level, computes force constants + from either displaced energies or displaced gradients. For energy-based + calculations, both gradients and Hessian elements are evaluated using + central finite differences. For gradient-based calculations, the Hessian + is obtained directly from gradient differences and subsequently + symmetrized. + + Returns + ------- + None + + Notes + ----- + The computed force-constant matrix is stored in ``self.FC``. For + energy-based calculations, the corresponding finite-difference + gradient is stored in ``self.gradient``. + """ indices = self.indices disp = self.disp - if self.coord_type_b == "cartesian": + if self.coord_type == "cartesian": size = self.indices[-1][0] + 1 cart_disp = np.zeros(size) for i in range(len(cart_disp)): @@ -53,13 +178,12 @@ def run(self): else: denom_disp = self.disp.disp dim = self.p_array.shape[0] - self.gradient = np.zeros((dim)) self.FC = np.zeros((dim, dim)) if not self.deriv_level: + self.gradient = np.zeros((dim)) p_en_array = self.p_array m_en_array = self.m_array e_r = self.ref_en - print("Force Constant disp sizes:") for index in indices: i, j = index[0], index[1] e_pi, e_pj = p_en_array[i, i], p_en_array[j, j] @@ -85,21 +209,93 @@ def run(self): il = (cf[1], cf[0]) self.FC[il] = self.FC[cf] elif self.deriv_level == 1: - self.FC = (self.p_array - self.m_array) / (2 * denom_disp[0]) + self.FC = self.first_deriv(self.p_array, self.m_array, denom_disp[0]) + for i in range(len(self.FC) - 1): + for j in range(i + 1): + self.FC[i + 1, j] = (self.FC[i + 1, j] + self.FC[j, i + 1]) / 2 + self.FC[j, i + 1] = self.FC[i + 1, j] else: print("Higher order deriv_level computations aren't yet supported") raise RuntimeError # I might want to shift the above computation in the deriv_level == 1 down here, though it is only one line def first_deriv(self, e_p, e_m, disp): + """ + Compute a first derivative using a central finite difference. + + Parameters + ---------- + e_p : float or ndarray + Quantity evaluated at the positive displacement. + + e_m : float or ndarray + Quantity evaluated at the negative displacement. + + disp : float + Displacement magnitude. + + Returns + ------- + float or ndarray + First derivative estimate. + """ return (e_p - e_m) / (2 * disp) # Functions for computing the diagonal and off-diagonal force constants def diag_fc(self, e_p, e_m, e_r, disp): + """ + Compute a diagonal Hessian element by finite differences. + + Parameters + ---------- + e_p, e_m : float + Energies at positive and negative displacements. + + e_r : float + Reference energy. + + disp : float + Displacement magnitude. + + Returns + ------- + float + Diagonal force constant. + """ fc = (e_p - 2 * e_r + e_m) / (disp**2) return fc def off_diag_fc(self, e_pi_pj, e_pi, e_pj, e_mi, e_mj, e_mi_mj, e_r, disp1, disp2): + """ + Compute an off-diagonal Hessian element by finite differences. + + Parameters + ---------- + e_pi_pj : float + Energy at the simultaneous positive displacement of coordinates + ``i`` and ``j``. + + e_pi, e_pj : float + Energies at individual positive displacements. + + e_mi, e_mj : float + Energies at individual negative displacements. + + e_mi_mj : float + Energy at the simultaneous negative displacement of coordinates + ``i`` and ``j``. + + e_r : float + Reference energy. + + disp1, disp2 : float + Displacement magnitudes for coordinates ``i`` and ``j``. + + Returns + ------- + float + Off-diagonal force constant. + """ fc = (e_pi_pj - e_pi - e_pj + 2 * e_r - e_mi - e_mj + e_mi_mj) / ( 2 * disp1 * disp2 ) diff --git a/concordantmodes/g_matrix.py b/concordantmodes/g_matrix.py index f9d6c99..d7fae49 100644 --- a/concordantmodes/g_matrix.py +++ b/concordantmodes/g_matrix.py @@ -3,19 +3,126 @@ from numpy import linalg as LA -class GMatrix(object): +class GMatrix: """ - This class will be used to compute the G-matrix of the GF-Matrix Method. - It should be constructed to give the flexibility of transformation via - the L-matrix. + Construct the Wilson G-matrix for vibrational analysis. + + This class computes the kinetic energy matrix (G-matrix) used in the + Wilson GF method for determining molecular vibrational frequencies and + normal modes. The G-matrix is constructed from the Wilson B-matrix and + the inverse atomic mass matrix according to + + .. math:: + + G = B M^{-1} B^T + + where: + + - ``B`` is the Wilson B-matrix relating internal-coordinate + displacements to Cartesian displacements. + - ``M`` is the diagonal Cartesian mass matrix. + + The resulting G-matrix represents the kinetic energy metric in the + internal-coordinate basis and is subsequently combined with the + force-constant matrix in the generalized eigenvalue problem + + .. math:: + + GF L = L \\Lambda + + to obtain harmonic vibrational frequencies and normal coordinates. + + The class also supports projection into a reduced coordinate space, + such as delocalized internal coordinates or projected vibrational + coordinates used in the Concordant Modes Algorithm (CMA). + + Parameters + ---------- + zmat : Zmat + Molecular coordinate object containing atomic masses, atom labels, + and internal-coordinate definitions. + + s_vectors : SVectors + Object containing the Wilson B-matrix and related internal-coordinate + transformation information. + + options : object + User-defined options controlling coordinate-system selection and + projection behavior. + + proj : ndarray, optional + Projection matrix defining a reduced vibrational coordinate basis. + If provided, the G-matrix is transformed into the projected space. + + Attributes + ---------- + G : ndarray + Wilson G-matrix in either the full internal-coordinate basis or + the projected coordinate basis. + + zmat : Zmat + Molecular coordinate representation. + + s_vectors : SVectors + Internal-coordinate derivative object containing the Wilson + B-matrix. + + proj : ndarray + Projection matrix used to reduce the coordinate basis. + + Notes + ----- + Atomic masses are converted to inverse masses prior to constructing + the Cartesian mass matrix. Dummy atoms labeled ``"X"`` are assigned + zero inverse mass so that they do not contribute to the kinetic energy. + + Small numerical elements are removed in two stages: + + 1. Absolute thresholding to eliminate numerical noise. + 2. Relative thresholding based on the largest value in each row. + + This improves numerical stability for subsequent GF matrix + diagonalizations. """ - def __init__(self, zmat, s_vectors, options): + def __init__(self, zmat, s_vectors, options, proj=np.array([])): self.zmat = zmat self.s_vectors = s_vectors self.options = options + self.proj = proj def run(self): + """ + Compute the Wilson G-matrix. + + Constructs the inverse Cartesian mass matrix from the atomic masses, + evaluates + + .. math:: + + G = B M^{-1} B^T + + and optionally projects the result into a reduced coordinate basis. + + Small numerical elements are removed to improve numerical stability + during subsequent diagonalization of the GF matrix. + + Returns + ------- + None + + Notes + ----- + The computed matrix is stored in the instance attribute ``G``. + If a projection matrix is supplied and the coordinate system is not + a conventional Z-matrix representation, the projected matrix + + .. math:: + + G' = P^T G P + + is formed and stored instead. + """ B = self.s_vectors.B # Compute and temper G @@ -30,3 +137,16 @@ def run(self): u = np.diag(u) self.G = B.dot(u.dot(B.T)) self.G[np.abs(self.G) < tol] = 0 + + if len(self.proj) and self.options.coords != "ZMAT": + self.G = np.dot(self.proj.T, np.dot(self.G, self.proj)) + + del_tol = 1.0e-3 + for row in self.G: + abs_row = np.abs(row) + row[abs_row < np.max(abs_row) * del_tol] = 0 + self.G = self.G.T + for row in self.G: + abs_row = np.abs(row) + row[abs_row < np.max(abs_row) * del_tol] = 0 + self.G = self.G.T diff --git a/concordantmodes/g_read.py b/concordantmodes/g_read.py index ce0a263..81d6ddb 100644 --- a/concordantmodes/g_read.py +++ b/concordantmodes/g_read.py @@ -4,11 +4,64 @@ from numpy import linalg as LA -class GrRead(object): +class GrRead: """ - This class may be used to read in force constant information from the - standard 3-column formatted file. Idk what that format is called, but - it's pretty standard. + Read Cartesian gradient data from a formatted gradient file. + + This class extracts numerical gradient components from a text file and + reconstructs the Cartesian gradient vector corresponding to a molecular + geometry. The reader is intended for use in finite-difference Hessian + calculations, where gradients computed by an electronic structure + program are used to construct force-constant matrices. + + All floating-point values are extracted from the file using a regular + expression. The Cartesian gradient is assumed to be stored as the final + ``3N`` values in the file, where ``N`` is the number of atoms in the + molecular geometry. + + Parameters + ---------- + file_name : str + Name of the gradient file to be read. + + Attributes + ---------- + file_name : str + Input filename. + + grad_regex : re.Pattern + Regular expression used to extract floating-point values from the + gradient file. + + grad : list[str] + Raw numerical values extracted from the file. + + cart_grad : ndarray + Cartesian gradient vector containing the final ``3N`` gradient + components. + + Notes + ----- + The reader assumes that the gradient information appears at the end of + the file and that the total number of Cartesian components is equal to + + .. math:: + + 3N + + where ``N`` is the number of atoms. + + This utility is primarily used during gradient-based finite-difference + Hessian construction within the Concordant Modes Algorithm (CMA) + workflow. + + Examples + -------- + >>> reader = GrRead("fc_b.grad") + >>> reader.run(cartesians) + >>> grad = reader.cart_grad + >>> grad.shape + (3 * len(cartesians),) """ def __init__(self, file_name): @@ -16,6 +69,29 @@ def __init__(self, file_name): self.grad_regex = re.compile(r"(-?\d+\.\d+)") def run(self, cartesians): + """ + Read and extract the Cartesian gradient vector. + + Parameters + ---------- + cartesians : ndarray + Cartesian coordinates of the molecular geometry. The number of + coordinates is used to determine the expected length of the + Cartesian gradient vector. + + Returns + ------- + None + + Notes + ----- + All floating-point values are extracted from the file, and the final + ``3N`` values are interpreted as the Cartesian gradient components, + where ``N`` is the number of atoms. + + The resulting gradient vector is stored in the attribute + ``cart_grad``. + """ with open(self.file_name, "r") as file: grad = file.read() self.grad = re.findall(self.grad_regex, grad) diff --git a/concordantmodes/gf_method.py b/concordantmodes/gf_method.py index b415346..51d7eee 100644 --- a/concordantmodes/gf_method.py +++ b/concordantmodes/gf_method.py @@ -7,11 +7,146 @@ import copy -class GFMethod(object): +class GFMethod: """ - This class is a versatile method to compute GF frequencies. + Solve the Wilson GF vibrational eigenvalue problem. - TODO: Insert standard uncertainties of amu_elmass and HARTREE_WAVENUM + This class performs harmonic vibrational analysis using the Wilson + GF method. Given a force-constant matrix (F) and a kinetic energy + matrix (G), the generalized eigenvalue problem + + .. math:: + + GF L = L \\Lambda + + is transformed into an ordinary symmetric eigenvalue problem and + solved to obtain vibrational frequencies and normal coordinates. + + The procedure consists of: + + 1. Constructing the symmetric orthogonalizer + + .. math:: + + G^{1/2} + + 2. Forming the symmetrized force-constant matrix + + .. math:: + + F_O = G^{1/2} F G^{1/2} + + 3. Diagonalizing the transformed matrix + + .. math:: + + F_O L' = L' \\Lambda + + 4. Recovering the normal-coordinate transformation matrix + + .. math:: + + L = G^{1/2}L' + + 5. Converting the eigenvalues into harmonic vibrational frequencies. + + Optional molecular symmetry support allows the transformed force + matrix to be block diagonalized by irreducible representation, + preserving symmetry labels and reducing computational cost. + + Parameters + ---------- + G : ndarray + Wilson G-matrix describing the kinetic energy metric in the + chosen coordinate basis. + + F : ndarray + Harmonic force-constant matrix in the same coordinate basis + as the G-matrix. + + zmat : Zmat + Molecular coordinate object. + + ted : TED + Total Energy Distribution (TED) object used to analyze and + print vibrational mode compositions. + + options : object + User-defined options controlling symmetry handling, + numerical tolerances, and output formatting. + + symtext : object, optional + Symmetry information generated by the molecular symmetry + analysis routines. Required when symmetry block + diagonalization is enabled. + + cma : bool, optional + Indicates whether the calculation is being performed within + the Concordant Modes Algorithm workflow. + + sym_sort : list, optional + User-defined symmetry sorting information. + + Attributes + ---------- + G : ndarray + Wilson G-matrix. + + F : ndarray + Harmonic force-constant matrix. + + G_O : ndarray + Symmetric orthogonalizer, :math:`G^{1/2}`. + + F_O : ndarray + Symmetrized force matrix. + + eig_v : ndarray + Eigenvalues of the GF problem. + + freq : ndarray + Harmonic vibrational frequencies in cm⁻¹. + + L_p : ndarray + Eigenvectors of the symmetrized force matrix. + + L : ndarray + Normal-coordinate transformation matrix. + + irrep_labels : list + Symmetry labels assigned to vibrational modes when + symmetry analysis is enabled. + + irrep_degen : list + Degeneracies associated with the vibrational modes. + + HARTREE_WAVENUM : float + Conversion factor from Hartree atomic units to cm⁻¹. + + Notes + ----- + Frequencies are obtained from the square roots of the GF + eigenvalues: + + .. math:: + + \\omega_i = \\sqrt{\\lambda_i} + + Imaginary frequencies are reported as negative values by + converting the imaginary component of the eigenvalue square root. + + The resulting frequencies are converted from atomic units to + wavenumbers using + + .. math:: + + 1\\ \\text{Hartree} + = + 219474.6313708\\ \\text{cm}^{-1} + + After solving the GF problem, a Total Energy Distribution (TED) + analysis is performed to determine the contribution of each + internal coordinate to every normal mode. """ def __init__(self, G, F, zmat, ted, options, symtext=None, cma=None, sym_sort=[]): @@ -27,6 +162,32 @@ def __init__(self, G, F, zmat, ted, options, symtext=None, cma=None, sym_sort=[] self.HARTREE_WAVENUM = 219474.6313708 def run(self): + """ + Solve the GF vibrational eigenvalue problem. + + Constructs the orthogonalized force matrix, performs the + diagonalization, computes harmonic frequencies, generates the + normal-coordinate transformation matrix, and performs a Total + Energy Distribution (TED) analysis. + + Returns + ------- + None + + Notes + ----- + If molecular symmetry is enabled, the transformed force matrix is + block diagonalized by irreducible representation before + diagonalization. Otherwise, the full matrix is diagonalized + directly. + + Results are stored in the attributes: + + - ``freq`` : vibrational frequencies + - ``L`` : normal-coordinate transformation matrix + - ``L_p`` : orthogonalized eigenvector matrix + - ``eig_v`` : GF eigenvalues + """ # Construct the orthogonalizer self.G_O = fractional_matrix_power(self.G, 0.5) @@ -90,6 +251,25 @@ def run(self): self.ted.run(self.L, self.freq, self.symtext, rect_print=False) def block_GF(self): + """ + Diagonalize the orthogonalized force matrix by symmetry blocks. + + The transformed force matrix is partitioned into blocks + corresponding to irreducible representations of the molecular + point group. Each block is diagonalized independently and the + resulting eigenvectors are assembled into a block-diagonal + transformation matrix. + + Returns + ------- + None + + Notes + ----- + This procedure preserves symmetry labels and mode ordering, + which is particularly useful for Concordant Modes calculations + that rely on symmetry-adapted coordinate transformations. + """ self.eigvals, self.eigvecs = [], [] offset_h = 0 self.block_fo = [] @@ -117,4 +297,23 @@ def block_GF(self): # function to extract bd matrices from supermatrix def extract(self, A, offset, size): + """ + Extract a square submatrix from a larger matrix. + + Parameters + ---------- + A : ndarray + Input matrix. + + offset : int + Starting row and column index. + + size : int + Dimension of the block to extract. + + Returns + ------- + ndarray + Extracted square submatrix. + """ return A[offset : offset + size, offset : offset + size] diff --git a/concordantmodes/int2cart.py b/concordantmodes/int2cart.py index cdf2383..e75be30 100644 --- a/concordantmodes/int2cart.py +++ b/concordantmodes/int2cart.py @@ -3,7 +3,7 @@ from . import masses -class Int2Cart(object): +class Int2Cart: def __init__(self, zmat, var_dict=None): self.zmat = zmat self.var_dict = var_dict diff --git a/concordantmodes/molden_writer.py b/concordantmodes/molden_writer.py index 946fdcd..f01890c 100644 --- a/concordantmodes/molden_writer.py +++ b/concordantmodes/molden_writer.py @@ -7,7 +7,73 @@ import shutil -class MoldenWriter(object): +class MoldenWriter: + """ + Write vibrational frequencies and normal-mode displacement vectors to + a Molden-format file. + + This class generates a Molden-compatible output file from the results + of a Concordant Modes calculation. The output contains the equilibrium + geometry, harmonic vibrational frequencies, and Cartesian displacement + vectors corresponding to the computed normal modes, allowing the modes + to be visualized and animated in Molden and other compatible molecular + visualization programs. + + Parameters + ---------- + zmat : Zmat + Molecular coordinate object containing atomic symbols and + Cartesian coordinates. + + disps : TransfDisp + Displacement object containing the normal-coordinate + displacement geometries. The diagonal elements of + ``p_disp`` are interpreted as positive displacements along + individual normal modes. + + freq : ndarray + Harmonic vibrational frequencies in cm⁻¹. + + Attributes + ---------- + zmat : Zmat + Molecular coordinate representation. + + disps : TransfDisp + Normal-mode displacement information. + + freq : ndarray + Harmonic vibrational frequencies. + + Notes + ----- + The generated Molden file contains the following sections: + + - ``[ATOMS]``: Atomic labels, atomic numbers, and Cartesian coordinates. + - ``[FREQ]``: Harmonic vibrational frequencies. + - ``[FR-COORD]``: Reference equilibrium geometry. + - ``[FR-NORM-COORD]``: Cartesian representations of the normal modes. + + The normal-mode vectors are obtained from the geometries generated by + displacements along individual normal coordinates. For mode ``i`` the + Cartesian displacement vector is computed as + + .. math:: + + \Delta \mathbf{r}_i = + \mathbf{r}_i^{(+)} - \mathbf{r}_0 + + where :math:`\mathbf{r}_0` is the reference geometry and + :math:`\mathbf{r}_i^{(+)}` is the geometry displaced along normal + coordinate :math:`Q_i`. + + These vectors therefore represent the Cartesian motion associated with + the vibrational normal modes obtained from the Concordant Modes + Algorithm. The vectors are multiplied by a constant scale factor prior + to writing the Molden file to improve visualization of vibrational + motion. + """ + def __init__( self, zmat, @@ -19,6 +85,30 @@ def __init__( self.freq = freq def run(self): + """ + Generate a Molden-format vibrational analysis file. + + Constructs the Molden sections describing the molecular geometry, + harmonic frequencies, reference coordinates, and normal-mode + displacement vectors, then writes the assembled output to a file + named ``MOLDEN``. + + Returns + ------- + None + + Notes + ----- + Atomic numbers are obtained using the QCEngine/QCElemental + periodic table utilities. + + Normal-mode displacement vectors are derived from the positive + displacement geometries stored in ``disps.p_disp`` and scaled + for improved visualization in Molden-compatible viewers. + + The generated file can be opened directly in Molden to animate + vibrational normal modes and inspect harmonic frequencies. + """ geometry = "" atom_num = np.array([]) for i in self.zmat.atom_list: diff --git a/concordantmodes/off_diag.py b/concordantmodes/off_diag.py index e7db137..60c3130 100644 --- a/concordantmodes/off_diag.py +++ b/concordantmodes/off_diag.py @@ -4,7 +4,7 @@ from scipy.linalg import fractional_matrix_power -class OffDiagonal(object): +class OffDiagonal: """ This class handles all off-diagonal CMA methods. """ diff --git a/concordantmodes/options.py b/concordantmodes/options.py index 0c954d2..b0e5412 100644 --- a/concordantmodes/options.py +++ b/concordantmodes/options.py @@ -1,4 +1,4 @@ -class Options(object): +class Options: def __init__(self, **kwargs): self.autosalcs = kwargs.pop("autosalcs", True) self.benchmark_full = kwargs.pop("benchmark_full", False) @@ -50,9 +50,9 @@ def __init__(self, **kwargs): # self.output_name_a = kwargs.pop("output_name_a", "output") self.partner_functions = kwargs.pop("partner_functions", True) self.printout_rel_e = kwargs.pop("printout_rel_e", True) - self.program_a = kwargs.pop("program_a", "molpro@2010.1.67+mpi") - self.program_b = kwargs.pop("program_b", "molpro@2010.1.67+mpi") - self.program_c = kwargs.pop("program_c", "molpro@2010.1.67+mpi") + self.program_a = kwargs.pop("program_a", "molpro") + self.program_b = kwargs.pop("program_b", "molpro") + self.program_c = kwargs.pop("program_c", "molpro") self.proj_tol = kwargs.pop("proj_tol", 1.0e-14) self.queue = kwargs.pop("queue", "batch") self.reduced_disp = kwargs.pop("reduced_disp", False) @@ -60,6 +60,7 @@ def __init__(self, **kwargs): self.second_order = kwargs.pop("second_order", False) self.rmsd = kwargs.pop("rmsd", False) self.ted_check = kwargs.pop("ted_check", False) + self.scaled_disp = kwargs.pop("scaled_disp", False) self.spin = kwargs.pop("spin", 1) self.subgroup = kwargs.pop("subgroup", False) self.success_regex_a = kwargs.pop("success_regex_a", "") diff --git a/concordantmodes/reap.py b/concordantmodes/reap.py index a18c7a2..059a6a8 100644 --- a/concordantmodes/reap.py +++ b/concordantmodes/reap.py @@ -1,13 +1,140 @@ import numpy as np -import json import os -import shutil import re +from pathlib import Path +from numpy.linalg import pinv from concordantmodes.s_vectors import SVectors -class Reap(object): +class Reap: + """ + Collect energies and gradients from displaced CMA calculations. + + This class extracts electronic energies or Cartesian gradients from + quantum chemistry calculations performed on displaced geometries + generated by the Concordant Modes Algorithm (CMA). The harvested + quantities are assembled into arrays suitable for finite-difference + construction of gradients and force-constant matrices. + + The class supports both: + + * Energy-based finite differences (``deriv_level=0``) + - Reads electronic energies from displaced calculations. + - Constructs arrays of positive and negative displacement energies. + - Stores the reference energy. + + * Gradient-based finite differences (``deriv_level=1``) + - Reads Cartesian gradients from displaced calculations. + - Constructs arrays of positive and negative displacement gradients. + - Stores the reference gradient. + + Program-specific output formats are handled through user-supplied + regular expressions, allowing the class to operate with multiple + electronic structure packages. + + Parameters + ---------- + options : object + User-defined options containing regular expressions, program + settings, symmetry options, and coordinate conversion controls. + + num_deg_free : int + Number of vibrational degrees of freedom. + + indices : list[list[int]] + Coordinate index pairs defining the displaced coordinates + generated by the CMA displacement algorithm. + + symm_obj : Symmetry + Symmetry object used to map displacement indices and exploit + point-group symmetry. + + cma_level : {"A", "B", "C"} + CMA stage associated with the calculations being harvested. + + deriv_level : int, optional + Differentiation level. + + - ``0`` : Read energies. + - ``1`` : Read gradients. + + disp : TransfDisp, optional + Displacement object containing displaced geometries and + coordinate transformation information. + + proj : ndarray, optional + Projection matrix defining the vibrational coordinate basis. + + zmat : Zmat, optional + Molecular coordinate object required when gradient conversion + into projected internal coordinates is requested. + + Attributes + ---------- + p_en_array : ndarray + Energies evaluated at positive displacements. + + m_en_array : ndarray + Energies evaluated at negative displacements. + + ref_en : float + Reference energy at the equilibrium geometry. + + p_grad_array : ndarray + Gradients evaluated at positive displacements. + + m_grad_array : ndarray + Gradients evaluated at negative displacements. + + ref_grad : ndarray + Gradient at the reference geometry. + + fail_list : list + Directories corresponding to failed electronic structure jobs. + + energy_re : re.Pattern + Compiled regular expression used to locate energies. + + success_re : re.Pattern + Compiled regular expression used to verify successful job + completion. + + Notes + ----- + The expected directory structure is + + :: + + DispsX/ + 1/ + 2/ + 3/ + ... + + where directory ``1`` contains the reference geometry and + subsequent directories contain positive and negative displacement + calculations. + + For gradient-based calculations, gradients may optionally be + transformed from Cartesian coordinates into the projected internal + coordinate basis using + + .. math:: + + A = B^{+} P + + where :math:`B^{+}` is the pseudoinverse of the Wilson B-matrix and + :math:`P` is the projection matrix defining the vibrational + coordinate space. + """ + + REGEX_MAP = { + "A": ("gradient_regex_a", "energy_regex_a", "success_regex_a"), + "B": ("gradient_regex_b", "energy_regex_b", "success_regex_b"), + "C": ("gradient_regex_c", "energy_regex_c", "success_regex_c"), + } + def __init__( self, options, @@ -16,9 +143,8 @@ def __init__( symm_obj, cma_level, deriv_level=0, - disp_sym=None, disp=None, - ted=None, + proj=None, zmat=None, ): self.options = options @@ -28,242 +154,263 @@ def __init__( self.deriv_level = deriv_level self.cma_level = cma_level self.disp = disp - self.ted = ted + self.proj = proj if proj is not None else np.array([]) self.zmat = zmat - if cma_level == "B": - if self.deriv_level: - self.gradient_regex = self.options.gradient_regex_b - self.energy_regex = self.options.energy_regex_b - self.success_regex = self.options.success_regex_b - else: - self.energy_regex = self.options.energy_regex_b - self.success_regex = self.options.success_regex_b - elif cma_level == "C": - if self.deriv_level: - pass - else: - self.energy_regex = self.options.energy_regex_c - self.success_regex = self.options.success_regex_c - else: # cma_level = "A" - if self.deriv_level: - self.gradient_regex = self.options.gradient_regex_a - self.energy_regex = self.options.energy_regex_a - self.success_regex = self.options.success_regex_a - else: - self.energy_regex = self.options.energy_regex_a - self.success_regex = self.options.success_regex_a + self.fail_list = [] + self.ref_grad = [] + + self._init_regex() + + def _init_regex(self): + try: + grad_key, energy_key, success_key = self.REGEX_MAP[self.cma_level] + except KeyError: + raise RuntimeError("cma_level must be A, B, or C") + + self.energy_regex = getattr(self.options, energy_key) + self.success_regex = getattr(self.options, success_key) + + if self.deriv_level: + self.gradient_regex = getattr(self.options, grad_key) + + # Compile once + self.energy_re = re.compile(self.energy_regex) + self.success_re = re.compile(self.success_regex) + + if self.deriv_level: + self.grad_re1 = re.compile(self.gradient_regex[0]) + self.grad_re2 = re.compile(self.gradient_regex[1]) + + # ------------------------- + # Core runner + # ------------------------- def run(self): - # Define energy/gradient search regex - if not self.deriv_level: - energy_regex = re.compile(self.energy_regex) - success_regex = re.compile(self.success_regex) - self.energies = np.array([]) - else: - grad_regex1 = re.compile(self.gradient_regex[0]) - grad_regex2 = re.compile(self.gradient_regex[1]) + """ + Harvest energies or gradients from displacement calculations. + + Executes either the energy-based or gradient-based collection + workflow depending on the value of ``deriv_level``. + + Returns + ------- + None + + Notes + ----- + Results are stored in instance attributes and subsequently used + by the :class:`ForceConstant` class to construct finite-difference + gradients and Hessians. + """ + if self.deriv_level: + return self._run_gradients() + return self._run_energies() + + # ------------------------- + # Energy workflow + # ------------------------- + def _run_energies(self): + """ + Collect energies from displaced geometries. + + Reads the reference energy and all positive and negative + displacement energies, storing them in square matrices indexed + according to the CMA displacement coordinates. + + Returns + ------- + None + + Notes + ----- + The resulting arrays ``p_en_array`` and ``m_en_array`` are used + for finite-difference evaluation of gradients and force constants. + """ + print(os.getcwd()) + print("Checking energies...") + + ref_en = self._read_energy(1) + print(f"Reference energy: {ref_en:.10f}") + size = self.num_deg_free + p_en = np.zeros((size, size)) + m_en = np.zeros((size, size)) - # if self.options.second_order: - # print(self.indices) - # print(self.indices[-1]) - # raise RuntimeError - # size = self.indices[-1][0] + 1 + indices = self._resolve_indices() - self.fail_list = [] - if not self.deriv_level: - print( - "If something looks wrong with the final frequencies, check these energies!" - ) - print("(Job number 1 == Reference energy) :D") - print(os.getcwd()) - os.chdir("./" + str(1)) - with open("output.dat", "r") as file: - data = file.read() - print(f"Success regex {success_regex}") - if not re.search(success_regex, data): - print("Energy failed at " + str("ref")) - raise RuntimeError - os.chdir("..") - - ref_en = float(re.findall(energy_regex, data)[0]) - print("Reference energy: " + str(ref_en)) - - indices = self.indices - p_en_array = np.zeros((size, size)) - m_en_array = np.zeros((size, size)) - rel_en_p = np.zeros((size, size)) - rel_en_m = np.zeros((size, size)) - relative_energies = [] - absolute_energies = [[("ref", "ref"), "ref", ref_en, 1]] - - direc = 2 - - if self.symm_obj.symtext is not None and self.options.exploit_pm_symm: - if self.options.only_TSIR: - print("Reap only the TSIR displacements") - indices = self.symm_obj.indices_by_irrep[0] - else: - print("Reap displacements from all irreps") - indices = self.symm_obj.indices_by_irrep.flatten().reshape((-1, 2)) - - Sum = 1 - h = 0 - # print(p_en_array.shape) - # print(p_en_array) - # print(indices) - for index in indices: - i, j = index[0], index[1] - p_en_array[i, j] = energy = self.reap_energies( - direc, success_regex, energy_regex - ) - print("p_en") - print(energy) - rel = energy - ref_en - print( - "Relative plus " - + "{:4d}".format(direc) - + "{:4d}".format(i) - + " " - + "{:4d}".format(j) - + ": " - + "{: 10.9f}".format(rel) - ) - rel_en_p[i, j] = rel - relative_energies.append([(i, j), "plus", rel, direc]) - absolute_energies.append([(i, j), "plus", energy, direc]) - - if h != 0: - m_en_array[i, j] = energy = self.reap_energies( - direc, success_regex, energy_regex - ) - direc += 1 - else: - m_en_array[i, j] = energy = self.reap_energies( - direc + 1, success_regex, energy_regex - ) - print("m_en") - print(energy) - rel = energy - ref_en - print( - "Relative minus " - + "{:4d}".format(direc + 1) - + "{:4d}".format(i) - + " " - + "{:4d}".format(j) - + ": " - + "{: 10.9f}".format(rel) - ) - rel_en_m[i, j] = rel - relative_energies.append([(i, j), "minus", rel, direc + 1]) - absolute_energies.append([(i, j), "minus", energy, direc + 1]) - direc += 2 - - if ( - self.symm_obj.symtext is not None - and self.options.exploit_pm_symm - and not self.options.only_TSIR - and Sum > len(self.symm_obj.indices_by_irrep[0]) - ): - h += 1 - Sum += 1 - - self.p_en_array = p_en_array - self.m_en_array = m_en_array - self.ref_en = ref_en - print_en = absolute_energies - np.set_printoptions(precision=2, linewidth=120) - os.chdir("..") - - else: - indices = self.indices - p_grad_array = np.array([], dtype=object) - m_grad_array = np.array([], dtype=object) - Sum = 0 - print(indices) - for index in indices: - grad = self.reap_gradients( - 2 * index[0] + 1 - Sum, grad_regex1, grad_regex2 - ) - p_grad_array = np.append(p_grad_array, grad, axis=0) - grad = self.reap_gradients( - 2 * index[0] + 2 - Sum, grad_regex1, grad_regex2 - ) - m_grad_array = np.append(m_grad_array, grad, axis=0) - self.p_grad_array = p_grad_array.reshape((-1, len(grad))) - self.m_grad_array = m_grad_array.reshape((-1, len(grad))) - - zmat_bool = self.zmat is not None - ted_bool = self.ted is not None - disp_bool = self.disp is not None - - if self.options.conv_grad and zmat_bool and ted_bool and disp_bool: - grad_s_vec = SVectors( - self.zmat, - self.options, - ) - - for i in range(len(indices)): - grad_s_vec.run(disp.p_disp[i], False) - A_proj = np.dot(LA.pinv(grad_s_vec.B), self.ted.proj) - self.p_array_b[i] = np.dot(cart_p_array_b[i].T, A_proj) - - grad_s_vec.run(disp.m_disp[i], False) - A_proj = np.dot(LA.pinv(grad_s_vec.B), self.ted.proj) - self.m_array_b[i] = np.dot(cart_m_array_b[i].T, A_proj) - - os.chdir("..") - if len(self.fail_list): - print("Some jobs have failed:") - print(self.fail_list) - raise RuntimeError - - def reap_energies(self, direc, success_regex, energy_regex): - os.chdir("./" + str(direc)) - - with open("output.dat", "r") as file: - data = file.read() - - if not re.search(success_regex, data): - self.fail_list.append(direc) - energy = 0.0 - else: - energy = float(re.findall(energy_regex, data)[0]) - - os.chdir("..") - - return energy - - def reap_gradients(self, direc, grad_regex1, grad_regex2): - os.chdir("./" + str(direc)) - grad_array = [] - if self.options.program_b == "molpro": - with open("output.xml", "r") as file: - data = file.readlines() - else: - with open("output.dat", "r") as file: - data = file.readlines() - for i in range(len(data)): - grad1 = re.search(grad_regex1, data[i]) - if grad1: - beg_grad = i + 1 - break - for i in range(len(data) - beg_grad): - grad2 = re.search(grad_regex2, data[i + beg_grad]) - if grad2: - end_grad = i + beg_grad - break - label_xyz = r"(\s*.*(\s*-?\d+\.\d+){3})+" - for line in data[beg_grad:end_grad]: - if re.search(label_xyz, line): - temp = line.split()[-3:] - grad_array.append(temp) - grad_array = np.array(grad_array) - grad_array = grad_array.astype("float64") - grad_array = grad_array.flatten() - if not grad1: - print("Gradient failed at " + os.getcwd()) - raise RuntimeError - os.chdir("..") - - return grad_array + direc = 2 + results = [] + + for i, j in indices: + e_plus = self._read_energy(direc) + p_en[i, j] = e_plus + results.append(((i, j), "plus", e_plus, direc)) + + e_minus = self._read_energy(direc + 1) + m_en[i, j] = e_minus + results.append(((i, j), "minus", e_minus, direc + 1)) + + print(f"[{i},{j}] Δ+ = {e_plus - ref_en:.9f}") + print(f"[{i},{j}] Δ- = {e_minus - ref_en:.9f}") + + direc += 2 + + self.p_en_array = p_en + self.m_en_array = m_en + self.ref_en = ref_en + + self._check_failures() + + # ------------------------- + # Gradient workflow + # ------------------------- + def _run_gradients(self): + """ + Collect gradients from displaced geometries. + + Reads Cartesian gradients from the reference and displaced + calculations and assembles them into arrays suitable for + gradient-based finite-difference Hessian construction. + + Returns + ------- + None + + Notes + ----- + If enabled through the options object, gradients may be converted + into the projected internal-coordinate basis before storage. + """ + p_list = [] + m_list = [] + + ref_grad = self._read_gradient(1) + + for idx in self.indices: + i = idx[0] + + p_grad = self._read_gradient(2 * i + 2) + m_grad = self._read_gradient(2 * i + 3) + + p_list.append(p_grad) + m_list.append(m_grad) + + self.ref_grad = np.array(ref_grad) + self.p_grad_array = np.array(p_list) + self.m_grad_array = np.array(m_list) + + if not self.options.cart_fc_b: + self._maybe_convert_gradients() + + self._check_failures() + + # ------------------------- + # Helpers + # ------------------------- + def _resolve_indices(self): + if self.symm_obj.symtext and self.options.exploit_pm_symm: + if self.options.only_TSIR: + return self.symm_obj.indices_by_irrep[0] + return self.symm_obj.indices_by_irrep.reshape((-1, 2)) + return self.indices + + def _read_energy(self, direc): + + direc = Path(str(direc)) + + data = self._read_file(direc, "output.dat") + + if not self.success_re.search(data): + self.fail_list.append(str(direc)) + return 0.0 + + return float(self.energy_re.findall(data)[0]) + + def _read_gradient(self, direc): + fname = "output.xml" if self.options.program_b == "molpro" else "output.dat" + lines = self._read_file(direc, fname, lines=True) + + start = next(i for i, l in enumerate(lines) if self.grad_re1.search(l)) + 1 + end = next( + i for i, l in enumerate(lines[start:], start) if self.grad_re2.search(l) + ) + + grad = [ + line.split()[-3:] + for line in lines[start:end] + if re.search(r"(\s*-?\d+\.\d+){3}", line) + ] + + return np.array(grad, dtype=float).flatten() + + def _read_file(self, direc, filename, lines=False): + path = Path(str(direc)) / filename + print(path) + with open(path, "r") as f: + return f.readlines() if lines else f.read() + + def _maybe_convert_gradients(self): + """ + Transform Cartesian gradients into the projected coordinate basis. + + For each displaced geometry, constructs the Wilson B-matrix and + computes the coordinate transformation + + .. math:: + + A = B^{+} P + + where :math:`B^{+}` is the Moore-Penrose pseudoinverse of the + B-matrix and :math:`P` is the projection matrix. + + Returns + ------- + None + + Notes + ----- + The transformed gradients replace the Cartesian gradients stored + in ``p_grad_array`` and ``m_grad_array``. + """ + if not ( + self.options.conv_grad + and self.zmat is not None + and len(self.proj) + and self.disp is not None + and not self.deriv_level + ): + return + + svec = SVectors(self.zmat, self.options) + + p_grad_buff = [] + m_grad_buff = [] + + for i in range(len(self.indices)): + svec.run(self.disp.p_disp[i], False) + A = pinv(svec.B) @ self.proj + p_grad_buff.append((self.p_grad_array[i].T @ A).T) + + svec.run(self.disp.m_disp[i], False) + A = pinv(svec.B) @ self.proj + m_grad_buff.append((self.m_grad_array[i].T @ A).T) + + self.p_grad_array = np.array(p_grad_buff) + self.m_grad_array = np.array(m_grad_buff) + + def _check_failures(self): + """ + Verify that all electronic structure calculations completed. + + Returns + ------- + None + + Raises + ------ + RuntimeError + If one or more displacement calculations failed according to + the user-defined success criterion. + """ + if self.fail_list: + raise RuntimeError(f"Failed jobs: {self.fail_list}") diff --git a/concordantmodes/rmsd.py b/concordantmodes/rmsd.py index ad36d8c..6e84c9a 100644 --- a/concordantmodes/rmsd.py +++ b/concordantmodes/rmsd.py @@ -4,7 +4,7 @@ import time -class RMSD(object): +class RMSD: def __init__(self): print("Nothing to init") diff --git a/concordantmodes/s_vectors.py b/concordantmodes/s_vectors.py index 1629489..c8e3abb 100644 --- a/concordantmodes/s_vectors.py +++ b/concordantmodes/s_vectors.py @@ -1,18 +1,101 @@ -import os -import shutil import numpy as np from concordantmodes.ted import TED from concordantmodes.int2cart import Int2Cart -from numpy.linalg import inv from numpy import linalg as LA -class SVectors(object): +class SVectors: """ - s-vectors: s_a^k = (B_ax^k,B_ay^k,B_az^k) - Where a refers to the atom #, and k refers to the internal coordinate. - So, all B-tensors can be contained within the s-vector - for each atom in the molecular system of interest. + Construct first- and second-order Wilson B-tensors (s-vectors) for a + molecular internal coordinate system. + + This class evaluates the derivatives of internal coordinates with respect + to Cartesian coordinates and assembles them into the Wilson B-matrix used + throughout the GF method, coordinate transformations, force-constant + conversions, and vibrational analyses. Supported internal coordinate types + include bond stretches, angle bends, torsions, out-of-plane bends, linear + bends, and specialized linear bending coordinates (LinX and LinY). + + For an internal coordinate q_k and atom a, the corresponding s-vector is + + s_a^k = (∂q_k/∂x_a, ∂q_k/∂y_a, ∂q_k/∂z_a) + + and each row of the Wilson B-matrix is formed by concatenating the + s-vectors for all atoms associated with a given internal coordinate. + + The class can optionally: + + * Generate a projection matrix defining a linearly independent internal + coordinate basis using singular value decomposition (SVD). + * Compute numerical second-order B-tensors (B2) through finite + differentiation of first-order B-matrices. + * Support projected internal coordinate spaces for delocalized or + non-redundant coordinate analyses. + + Parameters + ---------- + zmat : Zmat + Molecular coordinate definition containing atom labels, masses, + Cartesian coordinates, and internal coordinate index lists + (bonds, angles, torsions, out-of-plane bends, etc.). + options : Options + Program options controlling coordinate generation, projection, + displacement magnitudes, numerical tolerances, and second-order + derivative calculations. + + Attributes + ---------- + bond_indices : list + Bond stretch coordinate definitions. + angle_indices : list + Angle bend coordinate definitions. + torsion_indices : list + Dihedral angle coordinate definitions. + oop_indices : list + Out-of-plane bending coordinate definitions. + lin_indices : list + Linear bend coordinate definitions. + linx_indices : list + Linear bend X-component coordinate definitions. + liny_indices : list + Linear bend Y-component coordinate definitions. + B : ndarray + First-order Wilson B-matrix with dimensions + (N_internal, 3N_atoms). + B2 : ndarray, optional + Numerical second-order Wilson B-tensor with dimensions + (N_internal, 3N_atoms, 3N_atoms). + proj : ndarray + Projection matrix mapping redundant internal coordinates to a + linearly independent coordinate space. + s_2center_dict : dict + First-order s-vectors for two-center coordinates (bonds). + s_3center_dict : dict + First-order s-vectors for three-center coordinates (angles). + s_4center_dict : dict + First-order s-vectors for four-center coordinates (torsions, + out-of-plane bends, and linear bends). + s_ncenter_dict : dict + Reserved storage for generalized multi-center coordinates. + + Notes + ----- + The analytical expressions for bond, angle, torsional, out-of-plane, + and linear-coordinate derivatives are based on the Wilson GF formalism. + Out-of-plane derivative expressions follow the conventions described in + Wilson, Decius, and Cross, *Molecular Vibrations*, with the atom ordering + adapted to the coordinate conventions used by Concordant Modes. + + When `second_order=True`, numerical differentiation of displaced + B-matrices is used to generate the second-order tensor B2, which is + required for non-stationary force constant transformations and + higher-order coordinate transformations. + + References + ---------- + Wilson, E. B.; Decius, J. C.; Cross, P. C. + *Molecular Vibrations: The Theory of Infrared and Raman Vibrational + Spectra*; Dover Publications, 1980. """ def __init__(self, zmat, options): @@ -28,10 +111,9 @@ def __init__(self, zmat, options): self.linx_indices = zmat.linx_indices self.liny_indices = zmat.liny_indices self.options = options - # self.variable_dict = variable_dict self.zmat = zmat - def run(self, carts, B_proj, proj=None, second_order=False): + def run(self, carts, B_proj, proj=np.array([]), second_order=False): self.proj = proj # Initialize the cartesian coordinates self.carts = carts @@ -48,494 +130,158 @@ def run(self, carts, B_proj, proj=None, second_order=False): # LinY Bends # First, bonds. - if len(self.bond_indices) > 0: - for i in range(len(self.bond_indices)): - indies = self.bond_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - - self.s_2center_dict["B" + str(i + 1)] = self.carts.copy() - self.s_2center_dict["B" + str(i + 1)] = ( - 0 * self.s_2center_dict["B" + str(i + 1)] - ) - b_1 = self.compute_e(x1, x2) - b_2 = -b_1 - if Len: - for j in indies[0]: - self.s_2center_dict["B" + str(i + 1)][int(j) - 1] += b_1 / len( - indies[0] - ) - for j in indies[1]: - self.s_2center_dict["B" + str(i + 1)][int(j) - 1] += b_2 / len( - indies[1] - ) - else: - self.s_2center_dict["B" + str(i + 1)][int(indies[0]) - 1] = b_1 - self.s_2center_dict["B" + str(i + 1)][int(indies[1]) - 1] = b_2 + for i, groups in enumerate(self.bond_indices): + x1, x2 = self._get_points(groups) + + e = self.compute_e(x1, x2) + values = [e, -e] + + key = f"B{i+1}" + self._build_s_vector(key, groups, values, self.s_2center_dict) # Next, angles. - if len(self.angle_indices) > 0: - for i in range(len(self.angle_indices)): - indies = self.angle_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - L = len(indies[0]) - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= L - x2 = [0.0, 0.0, 0.0] - L = len(indies[1]) - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= L - x3 = [0.0, 0.0, 0.0] - L = len(indies[2]) - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= L - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - - self.s_3center_dict["A" + str(i + 1)] = self.carts.copy() - self.s_3center_dict["A" + str(i + 1)] = ( - 0 * self.s_3center_dict["A" + str(i + 1)] - ) - r_1 = self.compute_r(x1, x2) - r_2 = self.compute_r(x2, x3) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x3, x2) - phi = self.compute_phi(e_1, e_2) - a_1 = self.compute_BEND(e_1, e_2, phi, r_1) - a_3 = self.compute_BEND(e_2, e_1, phi, r_2) - a_2 = -a_1 - a_3 - - if Len: - for j in indies[0]: - self.s_3center_dict["A" + str(i + 1)][int(j) - 1] += a_1 / len( - indies[0] - ) - for j in indies[1]: - self.s_3center_dict["A" + str(i + 1)][int(j) - 1] += a_2 / len( - indies[1] - ) - for j in indies[2]: - self.s_3center_dict["A" + str(i + 1)][int(j) - 1] += a_3 / len( - indies[2] - ) - else: - self.s_3center_dict["A" + str(i + 1)][int(indies[0]) - 1] = a_1 - self.s_3center_dict["A" + str(i + 1)][int(indies[1]) - 1] = a_2 - self.s_3center_dict["A" + str(i + 1)][int(indies[2]) - 1] = a_3 + for i, groups in enumerate(self.angle_indices): + x1, x2, x3 = self._get_points(groups) + + r1 = self.compute_r(x1, x2) + r2 = self.compute_r(x2, x3) + e1 = self.compute_e(x1, x2) + e2 = self.compute_e(x3, x2) + + phi = self.compute_phi(e1, e2) + + a1 = self.compute_BEND(e1, e2, phi, r1) + a3 = self.compute_BEND(e2, e1, phi, r2) + a2 = -a1 - a3 + + values = [a1, a2, a3] + + key = f"A{i+1}" + self._build_s_vector(key, groups, values, self.s_3center_dict) # Next, torsions. - if len(self.torsion_indices) > 0: - for i in range(len(self.torsion_indices)): - indies = self.torsion_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - x3 = [0.0, 0.0, 0.0] - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= len(indies[2]) - x4 = [0.0, 0.0, 0.0] - for j in indies[3]: - x4 += self.carts[int(j) - 1] - x4 /= len(indies[3]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - x4 = self.carts[int(indies[3]) - 1] - - self.s_4center_dict["D" + str(i + 1)] = self.carts.copy() - self.s_4center_dict["D" + str(i + 1)] = ( - 0 * self.s_4center_dict["D" + str(i + 1)] - ) - - r_1 = self.compute_r(x1, x2) - r_2 = self.compute_r(x2, x3) - r_3 = self.compute_r(x3, x4) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x2, x3) - e_3 = self.compute_e(x3, x4) - phi_1 = self.compute_phi(e_1, -e_2) - phi_2 = self.compute_phi(e_2, -e_3) - - t_1 = self.compute_TORS1(e_1, -e_2, phi_1, r_1) - t_4 = self.compute_TORS1(-e_3, e_2, phi_2, r_3) - t_2 = self.compute_TORS2(e_1, -e_2, -e_3, phi_1, phi_2, r_1, r_2) - t_3 = -t_1 - t_2 - t_4 - - if Len: - for j in indies[0]: - self.s_4center_dict["D" + str(i + 1)][int(j) - 1] += t_1 / len( - indies[0] - ) - for j in indies[1]: - self.s_4center_dict["D" + str(i + 1)][int(j) - 1] += t_2 / len( - indies[1] - ) - for j in indies[2]: - self.s_4center_dict["D" + str(i + 1)][int(j) - 1] += t_3 / len( - indies[2] - ) - for j in indies[3]: - self.s_4center_dict["D" + str(i + 1)][int(j) - 1] += t_4 / len( - indies[3] - ) - else: - self.s_4center_dict["D" + str(i + 1)][int(indies[0]) - 1] = t_1 - self.s_4center_dict["D" + str(i + 1)][int(indies[1]) - 1] = t_2 - self.s_4center_dict["D" + str(i + 1)][int(indies[2]) - 1] = t_3 - self.s_4center_dict["D" + str(i + 1)][int(indies[3]) - 1] = t_4 + for i, groups in enumerate(self.torsion_indices): + x1, x2, x3, x4 = self._get_points(groups) + + r_1 = self.compute_r(x1, x2) + r_2 = self.compute_r(x2, x3) + r_3 = self.compute_r(x3, x4) + e_1 = self.compute_e(x1, x2) + e_2 = self.compute_e(x2, x3) + e_3 = self.compute_e(x3, x4) + phi_1 = self.compute_phi(e_1, -e_2) + phi_2 = self.compute_phi(e_2, -e_3) + + t1 = self.compute_TORS1(e_1, -e_2, phi_1, r_1) + t4 = self.compute_TORS1(-e_3, e_2, phi_2, r_3) + t2 = self.compute_TORS2(e_1, -e_2, -e_3, phi_1, phi_2, r_1, r_2) + t3 = -t1 - t2 - t4 + + values = [t1, t2, t3, t4] + + key = f"D{i+1}" + self._build_s_vector(key, groups, values, self.s_4center_dict) # Now, out of plane bending. - if len(self.oop_indices) > 0: - for i in range(len(self.oop_indices)): - indies = self.oop_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - x3 = [0.0, 0.0, 0.0] - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= len(indies[2]) - x4 = [0.0, 0.0, 0.0] - for j in indies[3]: - x4 += self.carts[int(j) - 1] - x4 /= len(indies[3]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - x4 = self.carts[int(indies[3]) - 1] - - self.s_4center_dict["O" + str(i + 1)] = self.carts.copy() - self.s_4center_dict["O" + str(i + 1)] = ( - 0 * self.s_4center_dict["O" + str(i + 1)] - ) - - r_1 = self.compute_r(x1, x2) - r_2 = self.compute_r(x3, x2) - r_3 = self.compute_r(x4, x2) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x3, x2) - e_3 = self.compute_e(x4, x2) - phi = self.compute_phi(e_2, e_3) - theta = self.calc_OOP(x1, x2, x3, x4) - - o_1 = self.compute_OOP1(e_1, e_2, e_3, r_1, theta, phi) - o_3 = self.compute_OOP2(e_1, e_2, e_3, r_2, theta, phi) - o_4 = self.compute_OOP2(-e_1, e_3, e_2, r_3, theta, phi) - o_2 = -o_1 - o_3 - o_4 - - if Len: - for j in indies[0]: - self.s_4center_dict["O" + str(i + 1)][int(j) - 1] += o_1 / len( - indies[0] - ) - for j in indies[1]: - self.s_4center_dict["O" + str(i + 1)][int(j) - 1] += o_2 / len( - indies[1] - ) - for j in indies[2]: - self.s_4center_dict["O" + str(i + 1)][int(j) - 1] += o_3 / len( - indies[2] - ) - for j in indies[3]: - self.s_4center_dict["O" + str(i + 1)][int(j) - 1] += o_4 / len( - indies[3] - ) - else: - self.s_4center_dict["O" + str(i + 1)][int(indies[0]) - 1] = o_1 - self.s_4center_dict["O" + str(i + 1)][int(indies[1]) - 1] = o_2 - self.s_4center_dict["O" + str(i + 1)][int(indies[2]) - 1] = o_3 - self.s_4center_dict["O" + str(i + 1)][int(indies[3]) - 1] = o_4 + for i, groups in enumerate(self.oop_indices): + x1, x2, x3, x4 = self._get_points(groups) + + r_1 = self.compute_r(x1, x2) + r_2 = self.compute_r(x3, x2) + r_3 = self.compute_r(x4, x2) + e_1 = self.compute_e(x1, x2) + e_2 = self.compute_e(x3, x2) + e_3 = self.compute_e(x4, x2) + phi = self.compute_phi(e_2, e_3) + theta = self.calc_OOP(x1, x2, x3, x4) + + o1 = self.compute_OOP1(e_1, e_2, e_3, r_1, theta, phi) + o3 = self.compute_OOP2(e_1, e_2, e_3, r_2, theta, phi) + o4 = self.compute_OOP2(-e_1, e_3, e_2, r_3, theta, phi) + o2 = -o1 - o3 - o4 + + values = [o1, o2, o3, o4] + + key = f"O{i+1}" + self._build_s_vector(key, groups, values, self.s_4center_dict) # Linear bending. - if len(self.lin_indices) > 0: - for i in range(len(self.lin_indices)): - indies = self.lin_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - x3 = [0.0, 0.0, 0.0] - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= len(indies[2]) - x4 = [0.0, 0.0, 0.0] - for j in indies[3]: - x4 += self.carts[int(j) - 1] - x4 /= len(indies[3]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - x4 = self.carts[int(indies[3]) - 1] - - self.s_4center_dict["L" + str(i + 1)] = self.carts.copy() - self.s_4center_dict["L" + str(i + 1)] = ( - 0 * self.s_4center_dict["L" + str(i + 1)] - ) - r_1 = self.compute_r(x1, x2) - r_2 = self.compute_r(x3, x2) - r_3 = self.compute_r(x4, x2) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x3, x2) - e_3 = self.compute_e(x4, x2) - theta = self.calc_Lin(x1, x2, x3, x4) - - l_1 = self.compute_LIN(e_1, e_2, e_3, r_1, theta) - l_3 = self.compute_LIN(e_2, e_3, e_1, r_2, theta) - l_2 = -l_1 - l_3 - - if Len: - for j in indies[0]: - self.s_4center_dict["L" + str(i + 1)][int(j) - 1] += l_1 / len( - indies[0] - ) - for j in indies[1]: - self.s_4center_dict["L" + str(i + 1)][int(j) - 1] += l_2 / len( - indies[1] - ) - for j in indies[2]: - self.s_4center_dict["L" + str(i + 1)][int(j) - 1] += l_3 / len( - indies[2] - ) - else: - self.s_4center_dict["L" + str(i + 1)][int(indies[0]) - 1] = l_1 - self.s_4center_dict["L" + str(i + 1)][int(indies[1]) - 1] = l_2 - self.s_4center_dict["L" + str(i + 1)][int(indies[2]) - 1] = l_3 + for i, groups in enumerate(self.lin_indices): + x1, x2, x3, x4 = self._get_points(groups) + + r_1 = self.compute_r(x1, x2) + r_2 = self.compute_r(x3, x2) + r_3 = self.compute_r(x4, x2) + e_1 = self.compute_e(x1, x2) + e_2 = self.compute_e(x3, x2) + e_3 = self.compute_e(x4, x2) + theta = self.calc_Lin(x1, x2, x3, x4) + + l1 = self.compute_LIN(e_1, e_2, e_3, r_1, theta) + l3 = self.compute_LIN(e_2, e_3, e_1, r_2, theta) + l2 = -l1 - l3 + + values = [l1, l2, l3] + + key = f"L{i+1}" + self._build_s_vector(key, groups, values, self.s_4center_dict) # LinX bending. - if len(self.linx_indices) > 0: - for i in range(len(self.linx_indices)): - indies = self.linx_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - x3 = [0.0, 0.0, 0.0] - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= len(indies[2]) - x4 = [0.0, 0.0, 0.0] - for j in indies[3]: - x4 += self.carts[int(j) - 1] - x4 /= len(indies[3]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - x4 = self.carts[int(indies[3]) - 1] - - self.s_4center_dict["Lx" + str(i + 1)] = self.carts.copy() - self.s_4center_dict["Lx" + str(i + 1)] = ( - 0 * self.s_4center_dict["Lx" + str(i + 1)] - ) - r_1 = self.compute_r(x2, x1) - r_2 = self.compute_r(x3, x2) - r_3 = self.compute_r(x4, x3) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x2, x3) - e_3 = self.compute_e(x3, x4) - ax = self.calc_alpha_x(e_1, e_2, e_3) - phi_1 = self.compute_phi(-e_1, e_2) - phi_2 = self.compute_phi(-e_2, e_3) - - lx_1 = self.compute_LINX1(e_1, e_2, e_3, r_1, phi_1, phi_2, ax) - lx_2 = self.compute_LINX2(e_1, e_2, e_3, r_1, r_2, phi_1, phi_2, ax) - lx_4 = self.compute_LINX4(e_1, e_2, e_3, r_3, phi_1, ax) - lx_3 = -lx_1 - lx_2 - lx_4 - - if Len: - for j in indies[0]: - self.s_4center_dict["Lx" + str(i + 1)][ - int(j) - 1 - ] += lx_1 / len(indies[0]) - for j in indies[1]: - self.s_4center_dict["Lx" + str(i + 1)][ - int(j) - 1 - ] += lx_2 / len(indies[1]) - for j in indies[2]: - self.s_4center_dict["Lx" + str(i + 1)][ - int(j) - 1 - ] += lx_3 / len(indies[2]) - for j in indies[3]: - self.s_4center_dict["Lx" + str(i + 1)][ - int(j) - 1 - ] += lx_4 / len(indies[3]) - else: - self.s_4center_dict["Lx" + str(i + 1)][int(indies[0]) - 1] = lx_1 - self.s_4center_dict["Lx" + str(i + 1)][int(indies[1]) - 1] = lx_2 - self.s_4center_dict["Lx" + str(i + 1)][int(indies[2]) - 1] = lx_3 - self.s_4center_dict["Lx" + str(i + 1)][int(indies[3]) - 1] = lx_4 + for i, groups in enumerate(self.linx_indices): + x1, x2, x3, x4 = self._get_points(groups) + + r_1 = self.compute_r(x2, x1) + r_2 = self.compute_r(x3, x2) + r_3 = self.compute_r(x4, x3) + e_1 = self.compute_e(x1, x2) + e_2 = self.compute_e(x2, x3) + e_3 = self.compute_e(x3, x4) + ax = self.calc_alpha_x(e_1, e_2, e_3) + phi_1 = self.compute_phi(-e_1, e_2) + phi_2 = self.compute_phi(-e_2, e_3) + + lx1 = self.compute_LINX1(e_1, e_2, e_3, r_1, phi_1, phi_2, ax) + lx2 = self.compute_LINX2(e_1, e_2, e_3, r_1, r_2, phi_1, phi_2, ax) + lx4 = self.compute_LINX4(e_1, e_2, e_3, r_3, phi_1, ax) + lx3 = -lx1 - lx2 - lx4 + + values = [lx1, lx2, lx3, lx4] + + key = f"Lx{i+1}" + self._build_s_vector(key, groups, values, self.s_4center_dict) # LinY bending. - if len(self.liny_indices) > 0: - for i in range(len(self.liny_indices)): - indies = self.liny_indices[i] - Len = len(np.array(indies[0]).shape) - if Len: - x1 = [0.0, 0.0, 0.0] - for j in indies[0]: - x1 += self.carts[int(j) - 1] - x1 /= len(indies[0]) - x2 = [0.0, 0.0, 0.0] - for j in indies[1]: - x2 += self.carts[int(j) - 1] - x2 /= len(indies[1]) - x3 = [0.0, 0.0, 0.0] - for j in indies[2]: - x3 += self.carts[int(j) - 1] - x3 /= len(indies[2]) - x4 = [0.0, 0.0, 0.0] - for j in indies[3]: - x4 += self.carts[int(j) - 1] - x4 /= len(indies[3]) - else: - x1 = self.carts[int(indies[0]) - 1] - x2 = self.carts[int(indies[1]) - 1] - x3 = self.carts[int(indies[2]) - 1] - x4 = self.carts[int(indies[3]) - 1] - - self.s_4center_dict["Ly" + str(i + 1)] = self.carts.copy() - self.s_4center_dict["Ly" + str(i + 1)] = ( - 0 * self.s_4center_dict["Ly" + str(i + 1)] - ) - r_1 = self.compute_r(x2, x1) - r_2 = self.compute_r(x3, x2) - r_3 = self.compute_r(x4, x3) - e_1 = self.compute_e(x1, x2) - e_2 = self.compute_e(x2, x3) - e_3 = self.compute_e(x3, x4) - ay = self.calc_alpha_y(e_1, e_2, e_3) - phi_1 = self.compute_phi(-e_1, e_2) - phi_2 = self.compute_phi(-e_2, e_3) - - ly_1 = self.compute_LINY1(e_1, e_2, e_3, r_1, phi_1, ay) - ly_2 = self.compute_LINY2(e_1, e_2, e_3, r_1, r_2, phi_1, ay) - ly_4 = self.compute_LINY4(e_1, e_2, e_3, r_2, r_3, phi_1, ay) - ly_3 = -ly_1 - ly_2 - ly_4 - - if Len: - for j in indies[0]: - self.s_4center_dict["Ly" + str(i + 1)][ - int(j) - 1 - ] += ly_1 / len(indies[0]) - for j in indies[1]: - self.s_4center_dict["Ly" + str(i + 1)][ - int(j) - 1 - ] += ly_2 / len(indies[1]) - for j in indies[2]: - self.s_4center_dict["Ly" + str(i + 1)][ - int(j) - 1 - ] += ly_3 / len(indies[2]) - for j in indies[3]: - self.s_4center_dict["Ly" + str(i + 1)][ - int(j) - 1 - ] += ly_4 / len(indies[3]) - else: - self.s_4center_dict["Ly" + str(i + 1)][int(indies[0]) - 1] = ly_1 - self.s_4center_dict["Ly" + str(i + 1)][int(indies[1]) - 1] = ly_2 - self.s_4center_dict["Ly" + str(i + 1)][int(indies[2]) - 1] = ly_3 - self.s_4center_dict["Ly" + str(i + 1)][int(indies[3]) - 1] = ly_4 - - # The last step will be to concatenate all of the s-vectors into a single B-tensor, in the order shown below. - self.B = np.array([self.s_2center_dict["B1"].flatten()]) - - # Append stretches - for i in range(len(self.bond_indices) - 1): - self.B = np.append( - self.B, - np.array([self.s_2center_dict["B" + str(i + 2)].flatten()]), - axis=0, - ) + for i, groups in enumerate(self.liny_indices): + x1, x2, x3, x4 = self._get_points(groups) - # Append bends - for i in range(len(self.angle_indices)): - self.B = np.append( - self.B, - np.array([self.s_3center_dict["A" + str(i + 1)].flatten()]), - axis=0, - ) - # Append torsions - for i in range(len(self.torsion_indices)): - self.B = np.append( - self.B, - np.array([self.s_4center_dict["D" + str(i + 1)].flatten()]), - axis=0, - ) - # Append oop bends - for i in range(len(self.oop_indices)): - self.B = np.append( - self.B, - np.array([self.s_4center_dict["O" + str(i + 1)].flatten()]), - axis=0, - ) - # Append lin bends - for i in range(len(self.lin_indices)): - self.B = np.append( - self.B, - np.array([self.s_4center_dict["L" + str(i + 1)].flatten()]), - axis=0, - ) - # Append linx bends - for i in range(len(self.linx_indices)): - self.B = np.append( - self.B, - np.array([self.s_4center_dict["Lx" + str(i + 1)].flatten()]), - axis=0, - ) - # Append liny bends - for i in range(len(self.liny_indices)): - self.B = np.append( - self.B, - np.array([self.s_4center_dict["Ly" + str(i + 1)].flatten()]), - axis=0, - ) + r_1 = self.compute_r(x2, x1) + r_2 = self.compute_r(x3, x2) + r_3 = self.compute_r(x4, x3) + e_1 = self.compute_e(x1, x2) + e_2 = self.compute_e(x2, x3) + e_3 = self.compute_e(x3, x4) + ay = self.calc_alpha_y(e_1, e_2, e_3) + phi_1 = self.compute_phi(-e_1, e_2) + phi_2 = self.compute_phi(-e_2, e_3) + + ly1 = self.compute_LINY1(e_1, e_2, e_3, r_1, phi_1, ay) + ly2 = self.compute_LINY2(e_1, e_2, e_3, r_1, r_2, phi_1, ay) + ly4 = self.compute_LINY4(e_1, e_2, e_3, r_2, r_3, phi_1, ay) + ly3 = -ly1 - ly2 - ly4 + + values = [ly1, ly2, ly3, ly4] + + key = f"Ly{i+1}" + self._build_s_vector(key, groups, values, self.s_4center_dict) + + # The last step will be to concatenate all of the s-vectors into a single B-tensor. + rows = [] + + rows += [v.flatten() for v in self.s_2center_dict.values()] + rows += [v.flatten() for v in self.s_3center_dict.values()] + rows += [v.flatten() for v in self.s_4center_dict.values()] + + self.B = np.array(rows) tol = 1e-4 # Now we acquire a linearly independant set of internal coordinates from the diagonalized @@ -579,11 +325,7 @@ def run(self, carts, B_proj, proj=None, second_order=False): self.B2[i, k, j + 1] = self.B2[i, j + 1, k] self.B2 = self.B2.astype(float) self.proj = Proj - # np.set_printoptions(precision=6, linewidth=240) - # print(self.B) self.B = B - # print(self.B) - # raise RuntimeError def compute_STRE(self, x1, x2): s = (x1 - x2) / self.compute_r(x1, x2) @@ -802,3 +544,28 @@ def second_order_B(self): B2 = (B2 + np.swapaxes(B2, 1, 2)) / 2 return B2 + + def _get_points(self, groups): + """Return averaged coordinates for grouped or single indices.""" + pts = [] + for g in groups: + if np.ndim(g): # calculates centroid. May need logic for COM. + coords = np.mean([self.carts[int(i) - 1] for i in g], axis=0) + else: + coords = self.carts[int(g) - 1] + pts.append(coords) + return pts + + def _distribute(self, arr, groups, values): + """Distribute vector contributions across atoms.""" + for group, val in zip(groups, values): + if np.ndim(group): + for i in group: + arr[int(i) - 1] += val / len(group) + else: + arr[int(group) - 1] = val + + def _build_s_vector(self, key, groups, values, store): + arr = np.zeros_like(self.carts) + self._distribute(arr, groups, values) + store[key] = arr diff --git a/concordantmodes/sapelo_template.py b/concordantmodes/sapelo_template.py index 4c8fb39..f380946 100644 --- a/concordantmodes/sapelo_template.py +++ b/concordantmodes/sapelo_template.py @@ -3,7 +3,7 @@ from numpy import linalg as LA -class SapeloTemplate(object): +class SapeloTemplate: """ This file just stores the sapelo optstep script """ @@ -27,6 +27,7 @@ def __init__(self, options, job_num, prog_name, prog): "prog": prog, "tc": str(job_num), "cline": self.progdict[prog_name], + "silly": "{}", } # This can be inserted back in if the sync keyword is sorted # $ -sync y @@ -49,7 +50,7 @@ def __init__(self, options, job_num, prog_name, prog): export NSLOTS=4 export THREADS=1 -module load intel/2023a +module load intel/2022a # to change scratch dir to use local machine scratch export SCRATCH_DIR=/scratch/$USER/tmp/$SLURM_JOB_ID @@ -58,7 +59,7 @@ def __init__(self, options, job_num, prog_name, prog): export TMPDIR=$SCRATCH_DIR -mpirun -n $NSLOTS apptainer exec /work/jttlab/containers/molpro_mpipr.sif molpro.exe input.dat --output $SLURM_SUBMIT_DIR/output.dat --nouse-logfile --directory $SCRATCH_DIR +mpirun -n $NSLOTS apptainer exec /work/jttlab/containers/molpro-2021-gapr.sif molpro.exe input.dat --output $SLURM_SUBMIT_DIR/output.dat --nouse-logfile --directory $SCRATCH_DIR rm $SCRATCH_DIR -r @@ -100,17 +101,58 @@ def __init__(self, options, job_num, prog_name, prog): # rm $PSI_SCRATCH -r elif self.prog_name == "orca": self.sapelo_template = """#!/bin/bash -#SBATCH --job-name=Concordant # Job name (testBowtie2) -#SBATCH --partition=batch # Partition name (batch, highmem_p, or gpu_p) -#SBATCH --ntasks=10 # 1 task (process) for below commands -#SBATCH --cpus-per-task=1 # CPU core count per task, by default 1 CPU core per task -#SBATCH --mem-per-cpu=10G +#SBATCH --job-name=default # Job name +#SBATCH --partition=batch # Partition (queue) name +#SBATCH --constraint="EPYC|Intel" +#SBATCH --nodes=1 # Number of nodes +#SBATCH --ntasks=4 # Number of MPI ranks +#SBATCH --ntasks-per-node=4 # How many tasks on each node +#SBATCH --cpus-per-task=1 # Number of cores per MPI rank +#SBATCH --mem=10G # total Memory no more per processor #SBATCH --time={time_limit} -#SBATCH --output=%x_%j.out # Standard output log, e.g., testBowtie2_12345.out +#SBATCH --output="%x.%j".out # Standard output log +#SBATCH --error="%x.%j".err # Standard error log +#SBATCH --mail-user="%u"@uga.edu +#SBTACH --mail-type=END,FAIL + +cd $SLURM_SUBMIT_DIR +export NSLOTS=4 +export THREADS=1 # can try adjusting but I don't think orca uses much omp threading if any + +scratch_dir=/scratch/$USER/tmp/$SLURM_JOB_ID +mkdir -p $scratch_dir + +module=ORCA/6.1.0-OpenMPI-4.1.8-GCC-13.3.0-avx2 +module load $module +export OMP_NUM_THREADS=$THREADS + +# Set other variables +base=`basename input.dat .dat` + +# Copy Job/Executable Data +cp $SLURM_SUBMIT_DIR/input.dat $scratch_dir/input.dat +if [ -e $base.xyz ]; then cp $base.xyz $scratch_dir/guess.xyz ; fi +if [ -e $base.gbw ]; then cp $base.gbw $scratch_dir/guess.gbw ; fi +if [ -e $base.hess ]; then cp $base.hess $scratch_dir/guess.hess ; fi +if [ -e product.xyz ]; then cp product.xyz $scratch_dir/product.xyz ; fi +if [ -e ts_guess.xyz ]; then cp ts_guess.xyz $scratch_dir/ts_guess.xyz ; fi + +echo " Running orca on `hostname`" +echo " Running calculation..." + +cd $scratch_dir +/apps/eb/$module/bin/orca input.dat >& $SLURM_SUBMIT_DIR/output.dat || exit 1 + +echo " Saving data and cleaning up..." +# delete any temporary files that my be hanging around. +rm -f *.tmp* +find . -type f -size +50M -exec rm -f {silly} \; +tar --exclude='*tmp*' --transform "s,^,Job_Data_$SLURM_JOB_ID/," -vzcf $SLURM_SUBMIT_DIR/Job_Data_$SLURM_JOB_ID.tar.gz * + +echo " Job complete on `hostname`." + +rm $scratch_dir -r -ml OpenMPI/4.1.4-GCC-11.3.0 -ml ORCA/6.1.0-OpenMPI-4.1.8-GCC-13.3.0-avx2 -/apps/eb/ORCA/6.1.0-OpenMPI-4.1.8-GCC-13.3.0-avx2/bin/orca input.dat > output.dat """ def run(self): diff --git a/concordantmodes/sisyphus_template.py b/concordantmodes/sisyphus_template.py new file mode 100644 index 0000000..0817f88 --- /dev/null +++ b/concordantmodes/sisyphus_template.py @@ -0,0 +1,143 @@ +import numpy as np +from numpy.linalg import inv +from numpy import linalg as LA + + +class SisyphusTemplate: + """ + This file just stores the sapelo optstep script + """ + + def __init__(self, options, job_num, prog_name, prog): + self.prog_name = prog_name + self.progdict = { + "molpro": "molpro -n $NSLOTS --nouse-logfile --no-xml-output -o \ + output.dat input.dat", + "psi4": "psi4 -n $NSLOTS", + "cfour": prog + "+vectorization", + "xtb": "input.coord", + "gxtb": "input.coord", + } + self.odict = { + "q": options.queue, + "nslots": options.nslots, + "jarray": "1-{}".format(job_num), + "prog_name": prog_name, + "prog": prog, + "tc": str(job_num), + "cline": self.progdict[prog_name], + } + # This can be inserted back in if the sync keyword is sorted + # $ -sync y + if self.prog_name == "psi4": + self.sisyphus_template = """#!/bin/bash +#SBATCH --job-name=Concordant # Job name +#SBATCH --partition=batch # Partition (queue) name +#SBATCH --cpus-per-task=4 # Number of cores per MPI rank +#SBATCH --mem=40GB # Memory per processor +#SBATCH --time=72:00:00 +#SBATCH --output=output.dat +#SBATCH --hint=compute_bound + +source /nfs/cluster_config/modules.sh +module load psi4/nightly + +eval "$(conda shell.bash hook)" +# Load their personal Conda environment +conda activate p4dev +which psi4 +#capture the submission directory +SUBMIT_DIR="$SLURM_SUBMIT_DIR" + +#scratch directory +SCRATCH_DIR="/scratch/$USER/$SLURM_JOB_ID" + +srun --cpu-bind=verbose psi4 -i input.dat -o output.dat + + +echo "Job \$SLURM_JOB_ID running on \$HOSTNAME" +echo "CPU affinity:" +taskset -cp $$ +""" + + elif self.prog_name == "molpro": + self.sisyphus_template = """#!/bin/bash +#SBATCH --job-name=STEP--00-1 # Job name +#SBATCH --partition=batch # Partition (queue) name +#SBATCH --ntasks-per-node=8 # How many tasks on each node +#SBATCH --cpus-per-task=1 # Number of cores per MPI rank +#SBATCH --mem=85GB # Total memory +#SBATCH --time=10:00:00 +#SBATCH --output="%x.%j".out # Standard output log +#SBATCH --error="%x.%j".err # Standard error log + +# Load Molpro +module purge +module load molpro + +#Run Molpro + +export SCRATCH_DIR=/scratch/$USER/tmp/$SLURM_JOB_ID +mkdir -p $SCRATCH_DIR + +# Run Molpro +molpro -n $SLURM_NTASKS input.dat --output output.dat --directory $SCRATCH_DIR + +#ignored line --do not remove +""" + + elif self.prog_name == "gxtb": + self.sisyphus_template = """#!/bin/bash +#SBATCH --job-name=Concordant # Job name +#SBATCH --partition=batch # Partition (queue) name +#SBATCH --cpus-per-task=4 # Number of cores per MPI rank +#SBATCH --mem=40GB # Memory per processor +#SBATCH --time=72:00:00 +#SBATCH --output=output.dat +#SBATCH --hint=compute_bound + +#capture the submission directory +SUBMIT_DIR="$SLURM_SUBMIT_DIR" + +#scratch directory +SCRATCH_DIR="/scratch/$USER/$SLURM_JOB_ID" + +srun --cpu-bind=verbose gxtb -c input.coord + + +echo "Job \$SLURM_JOB_ID running on \$HOSTNAME" +echo "CPU affinity:" +taskset -cp $$ +""" + elif self.prog_name == "xtb": + self.sisyphus_template = """#!/bin/bash +#SBATCH --job-name=Concordant # Job name +#SBATCH --partition=batch # Partition (queue) name +#SBATCH --cpus-per-task=4 # Number of cores per MPI rank +#SBATCH --mem=40GB # Memory per processor +#SBATCH --time=72:00:00 +#SBATCH --output=output.dat +#SBATCH --hint=compute_bound + +source /nfs/cluster_config/modules.sh +module load xtb/6.7.1 + +#eval "$(conda shell.bash hook)" +## Load their personal Conda environment +#conda activate xtbenv +#capture the submission directory +SUBMIT_DIR="$SLURM_SUBMIT_DIR" + +#scratch directory +SCRATCH_DIR="/scratch/$USER/$SLURM_JOB_ID" + +srun --cpu-bind=verbose {xtb_variant} {cline} {xtb_addtl_opts} + + +echo "Job \$SLURM_JOB_ID running on \$HOSTNAME" +echo "CPU affinity:" +taskset -cp $$ +""" + + def run(self): + return self.sisyphus_template.format(**self.odict) diff --git a/concordantmodes/submit.py b/concordantmodes/submit.py index c649523..13ef197 100644 --- a/concordantmodes/submit.py +++ b/concordantmodes/submit.py @@ -8,9 +8,132 @@ from concordantmodes.vulcan_template import VulcanTemplate from concordantmodes.sapelo_template import SapeloTemplate +from concordantmodes.sisyphus_template import SisyphusTemplate -class Submit(object): +class Submit: + """ + Submit and monitor batches of Concordant Modes displacement calculations on + supported computing clusters. + + This class automates the execution of electronic structure calculations + associated with generated displacement geometries by creating appropriate + scheduler submission scripts, launching jobs, and optionally monitoring + their completion. Supported execution environments include Sun Grid Engine + (SGE), Slurm, and user-defined custom submission workflows. + + Depending on the selected cluster type, the class will: + + * Generate scheduler-specific submission scripts using predefined + templates. + * Submit displacement calculations as either job arrays (SGE) or + individual jobs (Slurm). + * Monitor submitted jobs until completion when supported. + * Support externally managed submission workflows through a user-supplied + command string. + * Return control to the CMA workflow once all displacement calculations + have completed successfully. + + Parameters + ---------- + options : Options + Program options object containing cluster configuration, scheduler + settings, submission commands, and resource specifications. + cma_level : str + Concordant Modes displacement level being processed ("A", "B", or + "C"). Used to locate the corresponding displacement directory. + rootdir : str + Root working directory containing the displacement calculation + directories. + prog_name : str + Name of the electronic structure package being executed (e.g. + Molpro, Psi4, CFOUR, ORCA). + prog : str + Program execution command or executable path used by the submission + templates. + + Attributes + ---------- + cma_level : str + Current CMA displacement level. + options : Options + Runtime options controlling submission behavior. + prog : str + Electronic structure program executable or launch command. + prog_name : str + Electronic structure package name. + rootdir : str + Root directory containing displacement calculations. + + Notes + ----- + The class assumes that displacement directories have already been created + and numbered sequentially beginning with directory ``1``. + + Scheduler-specific behavior: + + **SGE** + Generates a job-array submission script using ``VulcanTemplate``, + submits it with ``qsub``, and monitors completion through + ``qacct`` records. + + **Slurm** + Generates a submission script using ``SapeloTemplate``, copies the + script into each displacement directory, submits jobs using + ``sbatch``, and polls job status through ``sacct``. + + **Custom** + Executes a user-defined submission command specified by + ``options.custom_submit_str``. Because job completion cannot be + monitored generically, execution terminates after submission and + requires the user to rerun CMA after the jobs finish. + + Raises + ------ + RuntimeError + If an unsupported cluster type is specified. + RuntimeError + If the custom submission option is selected without a valid + submission command. + RuntimeError + After custom job submission to intentionally terminate execution + until external calculations have completed. + + See Also + -------- + VulcanTemplate + Generates SGE submission scripts. + SapeloTemplate + Generates Slurm submission scripts. + + Examples + -------- + Submit displacement calculations on a Slurm cluster:: + + submitter = Submit( + options, + cma_level="B", + rootdir=workdir, + prog_name="molpro", + prog="molpro" + ) + submitter.run() + + Submit calculations using a custom scheduler wrapper:: + + options.cluster = "custom" + options.custom_submit_str = "my_submit_command" + + submitter = Submit( + options, + cma_level="A", + rootdir=workdir, + prog_name="psi4", + prog="psi4" + ) + submitter.run() + """ + def __init__(self, options, cma_level, rootdir, prog_name, prog): self.cma_level = cma_level self.options = options @@ -24,6 +147,8 @@ def run(self): for i in os.listdir(self.rootdir + "/Disps" + self.cma_level): disp_list.append(i) + os.chdir(self.rootdir + "/Disps" + self.cma_level) + # TODO move Vulcan and Sapelo templates to more general sge and slurm templates. if self.options.cluster.lower() == "sge": v_template = VulcanTemplate( @@ -61,6 +186,9 @@ def run(self): s_template = SapeloTemplate( self.options, len(disp_list), self.prog_name, self.prog ) + # s_template = SisyphusTemplate( + # self.options, len(disp_list), self.prog_name, self.prog + # ) out = s_template.run() with open("sub_script.sh", "w") as file: @@ -126,7 +254,7 @@ def run(self): ) os.chdir("..") time.sleep(3) - os.chdir("..") + # os.chdir("..") print( "Jobs have been submitted. You will need to come back when they finish and run CMA again with relevent gen_disps and calc keywords set to false." @@ -135,6 +263,7 @@ def run(self): else: print( - "Only Vulcan, Sapelo, or Custom cluster options are available, select one of those!" + "Only SGE, Slurm, or Custom cluster options are available, select one of those!" ) raise RuntimeError + os.chdir("../") diff --git a/concordantmodes/sym_helper.py b/concordantmodes/sym_helper.py index c77e6b6..82820dd 100644 --- a/concordantmodes/sym_helper.py +++ b/concordantmodes/sym_helper.py @@ -4,7 +4,7 @@ from concordantmodes.zmat import Zmat -class SymHelper(object): +class SymHelper: """ This class handles sym_sort generation and helper functions for natural internal coordinates diff --git a/concordantmodes/symmetry.py b/concordantmodes/symmetry.py index e4450f6..904af47 100644 --- a/concordantmodes/symmetry.py +++ b/concordantmodes/symmetry.py @@ -16,7 +16,94 @@ raise ImportError("MolSym library not detected in your environment") -class Symmetry(object): +class Symmetry: + """ + Construct and manipulate symmetry-adapted linear combinations (SALCs) + of internal coordinates using the MolSym library. + + This class provides the interface between the Concordant Modes internal + coordinate representation and MolSym symmetry analysis. It generates + molecular symmetry information, constructs symmetry-adapted internal + coordinates, removes redundant coordinates within each irreducible + representation, and builds projection matrices suitable for symmetry- + blocked vibrational analyses. + + The resulting symmetry information may be used to: + + - Generate SALCs of internal coordinates. + - Project redundant internal coordinates into a nonredundant basis. + - Construct symmetry-blocked projection matrices. + - Reorder G and F matrices into symmetry block form. + - Exploit degeneracies and irreducible representations during + vibrational frequency calculations and Concordant Modes Analysis (CMA). + + Parameters + ---------- + zmat : Zmat + Molecular internal-coordinate representation containing atom + labels, Cartesian coordinates, coordinate definitions, and + indexing information. + options : Options + Runtime options controlling symmetry handling, subgroup + selection, degeneracy treatment, projection methods, and + coordinate generation. + proj : ndarray + User-supplied projection matrix. Used when + ``options.man_proj`` is enabled to define custom symmetry + coordinates. + + Attributes + ---------- + zmat : Zmat + Molecular coordinate object. + options : Options + User-defined runtime options. + proj : ndarray + Projection matrix used for manual symmetry adaptation. + symtext : molsym.Symtext + Molecular symmetry information generated by MolSym. + salcs : ProjectionOp + Symmetry-adapted linear combinations of internal coordinates. + new_proj : ndarray + Nonredundant symmetry projection matrix. + salc_proj : ndarray + Combined SALC projection matrix. + proj_irreps : list[int] + Number of nonredundant coordinates retained in each irreducible + representation. + salcblocks : list[ndarray] + SALC blocks grouped by irreducible representation. + sym_sort : list[list[int]] + Ordering information used to group coordinates according to + symmetry species. + + Notes + ----- + The internal-coordinate symmetry workflow proceeds as follows: + + 1. Build a QCSchema-like molecular description from the Z-matrix. + 2. Generate molecular symmetry information using MolSym. + 3. Construct internal-coordinate SALCs. + 4. Remove redundant coordinates within each symmetry block using SVD. + 5. Assemble a symmetry-adapted projection matrix. + 6. Optionally reorder G and F matrices into symmetry block form. + + When ``options.man_proj`` is enabled, user-generated projection + matrices are passed directly to MolSym's projection operator instead + of generating SALCs automatically. + + Warnings + -------- + MolSym currently contains known issues involving linear internal + coordinates. The implementation warns users when linear bending + coordinates are present because SALC generation may be unreliable + for such coordinate types. + + References + ---------- + MolSym: + https://github.com/NASymmetry/MolSym + """ def __init__(self, zmat, options, proj): self.zmat = zmat @@ -130,7 +217,7 @@ def package_salcs(self): self.salcs.salc_sets.append(ir_salcs) else: ir_salcs = [self.salcs[i].coeffs for i in self.salcs.salcs_by_irrep[ir]] - ir_salcs = np.row_stack(ir_salcs) + ir_salcs = np.vstack(ir_salcs) self.salcs.salc_sets.append(ir_salcs) fxn_list.append( [1 for i in range(0, (len(ir_salcs) // self.symtext.irreps[ir].d))] @@ -197,7 +284,7 @@ def make_proj(self, s_vec): newsblock.append(s) self.sblock = copy.deepcopy(sblock) - self.salc_proj = np.row_stack(newsblock).T + self.salc_proj = np.vstack(newsblock).T def create_flat_sym_sort(self, sym_sort): print("sym_sort: ") @@ -316,15 +403,6 @@ def mode_symmetry_sort(self, TED, sym_sort, freqs): sym_freqs = copy.deepcopy(sym_modes) del_list = [] for i in range(len(sym_modes)): - # if len(sym_modes[i]) > 1: - # for j in range(len(sym_modes[i])): - # index = sym_modes[i][j] - # sym_freqs[i][j] = freqs[index].copy() - # sym_freqs[i].reverse() - # elif len(sym_modes[i]) == 1: - # del_list.append(i) - # else: - # pass for j in range(len(sym_modes[i])): index = sym_modes[i][j] sym_freqs[i][j] = freqs[index].copy() @@ -333,8 +411,6 @@ def mode_symmetry_sort(self, TED, sym_sort, freqs): if len(del_list): for i in del_list: print(freqs[sym_modes[i][0]]) - # for i in del_list: - # del sym_freqs[i] flat_sym_freqs = [x for xs in sym_freqs for x in xs] flat_sym_freqs = np.array(flat_sym_freqs) diff --git a/concordantmodes/ted.py b/concordantmodes/ted.py index 3c60c1e..8dd811b 100644 --- a/concordantmodes/ted.py +++ b/concordantmodes/ted.py @@ -1,14 +1,90 @@ import numpy as np -from numpy.linalg import inv from numpy import linalg as LA -from scipy.linalg import fractional_matrix_power -class TED(object): +class TED: """ - This class will be used to print out total energy distributions - from internal coordinate eigenvalue matrices. If redundant - coordinates are utilized, projection matrices will be stored here. + Compute and print the Total Energy Distribution (TED) for vibrational modes. + + The TED quantifies the contribution of each internal coordinate to each + normal mode and is expressed as a percentage. Contributions are computed + from the normal-mode eigenvector matrix and its pseudoinverse according + to the standard TED formalism: + + TED_ij = L_ij * (L^-1)_ji * 100 + + where ``L`` is the normal-mode transformation matrix. The resulting TED + matrix provides a chemically intuitive decomposition of vibrational + motions into stretches, bends, torsions, out-of-plane bends, and linear + bending coordinates. + + Parameters + ---------- + proj : ndarray + Projection matrix used to transform eigenvectors from a reduced or + symmetry-adapted coordinate basis back into the full internal + coordinate representation. + zmat : Zmat + Molecular coordinate object containing internal-coordinate + definitions and indexing information used for TED labeling. + options : Options + Runtime options controlling output formatting and coordinate + handling. + + Attributes + ---------- + proj : ndarray + Projection matrix relating projected coordinates to the full + internal-coordinate basis. + zmat : Zmat + Molecular coordinate representation. + options : Options + User-defined runtime options. + symtext : molsym.Symtext or None + Molecular symmetry information used to annotate TED tables with + irreducible representation labels. + TED : ndarray + Total Energy Distribution matrix expressed as percentages. + + Notes + ----- + The TED matrix is computed from the projected normal-mode eigenvectors + and their Moore-Penrose pseudoinverse. Each column corresponds to a + vibrational normal mode and each row corresponds to an internal + coordinate. + + When ``rect_print=True``, the eigenvectors are first transformed back + into the full internal-coordinate basis using the projection matrix + before TED analysis. This is typically used when redundant or + symmetry-adapted coordinates have been employed during the vibrational + calculation. + + Coordinate labels are automatically generated from the Z-matrix + coordinate definitions using the following conventions: + + - B : Bond stretches + - A : Bond angles + - D : Dihedral/torsional angles + - O : Out-of-plane bends + - L : Linear bends + - Lx : Linear bend x-components + - Ly : Linear bend y-components + + If molecular symmetry information is available, TED tables include + point-group and irreducible-representation assignments for each + vibrational mode. + + Examples + -------- + >>> ted = TED(proj, zmat, options) + >>> ted.run(L_matrix, frequencies) + + Print TED using symmetry-adapted normal modes: + + >>> ted.run(L_matrix, frequencies, symtext=symtext) + + The resulting table reports the percentage contribution of each + internal coordinate to every vibrational mode. """ def __init__(self, proj, zmat, options): @@ -18,198 +94,99 @@ def __init__(self, proj, zmat, options): def run(self, eigs, freq, symtext=None, rect_print=True): self.symtext = symtext - proj_eigs = eigs - if len(np.shape(self.proj)) > 2 and np.shape(self.proj)[0] > 1: - proj_buff = self.proj[0] - for i in range(len(self.proj) - 1): - proj_buff = np.append(proj_buff, self.proj[i + 1], axis=1) - if rect_print: - proj_eigs = np.dot(self.proj, proj_eigs) + + proj_eigs = np.dot(self.proj, eigs) if rect_print else eigs proj_eigs_inv = LA.pinv(proj_eigs) + print("The eigenvectors (check these for phase):") print(proj_eigs) + self.TED = np.multiply(proj_eigs, proj_eigs_inv.T) * 100 self.table_print(freq, self.TED, rect_print) def table_print(self, freq, TED, rect_print): - if len(freq) != len(TED[0]): + if len(freq) != TED.shape[1]: print( - "Something has gone terribly wrong. Your Total Energy \ - distribution should have the same number of columns as \ - frequencies, but it does not." + "Something has gone terribly wrong. Your Total Energy " + "distribution should have the same number of columns as " + "frequencies, but it does not." ) - print("TED columns: " + str(len(TED[0]))) + print("TED columns: " + str(TED.shape[1])) print("# frequencies: " + str(len(freq))) raise RuntimeError + div = 15 - a = len(freq) // div - for i in range(a): - print( - self.sub_table(freq[i * div : (i + 1) * div], TED, i, div, rect_print) - ) - if len(freq) % div: - print(self.sub_table(freq[a * div :], TED, a, div, rect_print)) + chunks = [freq[i : i + div] for i in range(0, len(freq), div)] + + for idx, chunk in enumerate(chunks): + print(self.sub_table(chunk, TED, idx, div, rect_print)) + + def _coord_label(self, i): + """Return formatted coordinate label prefix (exact same formatting).""" + z = self.zmat + + bounds = [ + (len(z.bond_indices), "B", " STRE: "), + (len(z.angle_indices), "A", " BEND: "), + (len(z.torsion_indices), "D", " TORS: "), + (len(z.oop_indices), "O", " OOP: "), + (len(z.lin_indices), "L", " LIN: "), + (len(z.linx_indices), "Lx", " LINX: "), + (len(z.liny_indices), "Ly", " LINY: "), + ] + + offset = 0 + for size, label, suffix in bounds: + if i < offset + size: + k = i - offset + return ( + "{:15s}".format(" ") + "{:4s}".format(label + str(k + 1)) + suffix + ) + offset += size + + return "{:15s}".format(" ") def sub_table(self, freq, TED, int_div, div, rect_print): - table_output = "{:>26s}".format("frequency #: ") n = len(freq) + + table_output = "{:>26s}".format("frequency #: ") for l in range(n): table_output += "{:8d}".format(l + int_div * div + 1) table_output += "\n" - table_output += "--------------------------" - for l in range(n): - table_output += "--------" - table_output += "\n" + + table_output += "--------------------------" + "--------" * n + "\n" + irrep_labels = [] if self.symtext is not None: table_output += "Point Group " + str(self.symtext.pg) + " " for ir, h in enumerate(self.symtext.salcblocks): - for s in range(h.shape[0]): - irrep_labels.append(str(self.symtext.irreps[ir].symbol)) + irrep_labels.extend([str(self.symtext.irreps[ir].symbol)] * h.shape[0]) table_output += "\n" + table_output += "{:>26s}".format("frequency: ") - for l in range(n): - table_output += " " + "{:7.1f}".format(freq[l]) + for val in freq: + table_output += " " + "{:7.1f}".format(val) table_output += "\n" + if self.symtext is not None: table_output += "{:>26s}".format("Irreps: ") for l in range(n): table_output += " " + "{:>7}".format(irrep_labels[l]) - table_output += "\n" - table_output += "--------------------------" - for l in range(n): - table_output += "--------" - table_output += "\n" + table_output += "\n" + + table_output += "--------------------------" + "--------" * n + "\n" + if rect_print: for i in range(len(TED)): - for j in range(len(freq)): - if i < len(self.zmat.bond_indices) and j == 0: - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("B" + str(i + 1)) - + " STRE: " - ) - elif ( - i < len(self.zmat.bond_indices) + len(self.zmat.angle_indices) - and j == 0 - ): - k = i - len(self.zmat.bond_indices) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("A" + str(k + 1)) - + " BEND: " - ) - elif ( - i - < len(self.zmat.bond_indices) - + len(self.zmat.angle_indices) - + len(self.zmat.torsion_indices) - and j == 0 - ): - k = ( - i - - len(self.zmat.bond_indices) - - len(self.zmat.angle_indices) - ) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("D" + str(k + 1)) - + " TORS: " - ) - elif ( - i - < len(self.zmat.bond_indices) - + len(self.zmat.angle_indices) - + len(self.zmat.torsion_indices) - + len(self.zmat.oop_indices) - and j == 0 - ): - k = ( - i - - len(self.zmat.bond_indices) - - len(self.zmat.angle_indices) - - len(self.zmat.torsion_indices) - ) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("O" + str(k + 1)) - + " OOP: " - ) - elif ( - i - < len(self.zmat.bond_indices) - + len(self.zmat.angle_indices) - + len(self.zmat.torsion_indices) - + len(self.zmat.oop_indices) - + len(self.zmat.lin_indices) - and j == 0 - ): - k = ( - i - - len(self.zmat.bond_indices) - - len(self.zmat.angle_indices) - - len(self.zmat.torsion_indices) - - len(self.zmat.oop_indices) - ) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("L" + str(k + 1)) - + " LIN: " - ) - elif ( - i - < len(self.zmat.bond_indices) - + len(self.zmat.angle_indices) - + len(self.zmat.torsion_indices) - + len(self.zmat.oop_indices) - + len(self.zmat.lin_indices) - + len(self.zmat.linx_indices) - and j == 0 - ): - k = ( - i - - len(self.zmat.bond_indices) - - len(self.zmat.angle_indices) - - len(self.zmat.torsion_indices) - - len(self.zmat.oop_indices) - - len(self.zmat.lin_indices) - ) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("Lx" + str(k + 1)) - + " LINX: " - ) - elif ( - i - < len(self.zmat.bond_indices) - + len(self.zmat.angle_indices) - + len(self.zmat.torsion_indices) - + len(self.zmat.oop_indices) - + len(self.zmat.lin_indices) - + len(self.zmat.linx_indices) - + len(self.zmat.liny_indices) - and j == 0 - ): - k = ( - i - - len(self.zmat.bond_indices) - - len(self.zmat.angle_indices) - - len(self.zmat.torsion_indices) - - len(self.zmat.oop_indices) - - len(self.zmat.lin_indices) - - len(self.zmat.linx_indices) - ) - table_output += ( - "{:15s}".format(" ") - + "{:4s}".format("Ly" + str(k + 1)) - + " LINY: " - ) + for j in range(n): + if j == 0: + table_output += self._coord_label(i) table_output += "{:8.1f}".format(TED[i][j + div * int_div]) table_output += "\n" else: for i in range(len(TED)): - for j in range(len(freq)): - if not j: + for j in range(n): + if j == 0: table_output += ( "{:>21s}".format("Mode #") + "{:>3s}".format(str(i + 1)) diff --git a/concordantmodes/tests/_test_symmetry.py b/concordantmodes/tests/_test_symmetry.py deleted file mode 100644 index 7b91ed5..0000000 --- a/concordantmodes/tests/_test_symmetry.py +++ /dev/null @@ -1,88 +0,0 @@ -import os -from concordantmodes.symmetry import Symmetry -from concordantmodes.options import Options -from concordantmodes.zmat import Zmat -from concordantmodes.s_vectors import SVectors - -# import molsym -import numpy as np - - -def test_symmetry(): - os.chdir("./ref_data/symmetry_test/") - test_subjects_del = ["c2v_water_zmat_del", "c3v_ammonia_zmat_del"] - test_subjects_custom = ["c2v_water_zmat_custom", "c3v_ammonia_zmat_custom"] - # delocalized first - test_subject_projs_del = [ - np.array( - [ - [0.56086644, 0.43061449, -0.70710678], - [0.56086644, 0.43061449, 0.70710678], - [-0.60898085, 0.79318492, 0.0], - ] - ), - np.array( - [ - [-0.14852746, 0.55791839, 0.3291588, -0.24149772], - [-0.14852746, 0.55791839, 0.3291588, -0.24149772], - [-0.14852746, 0.55791839, -0.6583176, 0.48299545], - [0.23900037, 0.06362601, -0.42112908, -0.5739944], - [0.23900037, 0.06362601, 0.21056454, 0.2869972], - [0.23900037, 0.06362601, 0.21056454, 0.2869972], - [0.50413466, 0.13420931, 0.23650561, 0.3223546], - [-0.50413466, -0.13420931, 0.1182528, 0.1611773], - [0.50413466, 0.13420931, -0.1182528, -0.1611773], - ] - ), - ] - - # start of customs - test_subject_projs_custom = [ - np.array( - [ - [0.56086644, 0.43061449, -0.70710678], - [0.56086644, 0.43061449, 0.70710678], - [-0.60898085, 0.79318492, 0.0], - ] - ), - np.array( - [ - [-0.14852746, 0.55791839, 0.3291588, -0.24149772], - [-0.14852746, 0.55791839, 0.3291588, -0.24149772], - [-0.14852746, 0.55791839, -0.6583176, 0.48299545], - [0.23900037, 0.06362601, -0.42112908, -0.5739944], - [0.23900037, 0.06362601, 0.21056454, 0.2869972], - [0.23900037, 0.06362601, 0.21056454, 0.2869972], - [0.50413466, 0.13420931, -0.1182528, -0.1611773], - [0.50413466, 0.13420931, -0.1182528, -0.1611773], - [0.50413466, 0.13420931, 0.23650561, 0.3223546], - ] - ), - ] - # iterate over delocalized and custom coordinates - for c, coord in enumerate(["Delocalized", "Custom"]): - if coord == "Delocalized": - test_subjects = test_subjects_del - test_subject_projs = test_subject_projs_del - else: - test_subjects = test_subjects_custom - test_subject_projs = test_subject_projs_custom - - # iterate over zmats to test symmetry projection on - for t, test in enumerate(test_subjects): - options = Options() - options.cart_insert = 9 - options.molsym_symmetry = True - options.coords = coord - options.second_order = False - zmat = Zmat(options) - zmat.run(test) - symm_obj = Symmetry(zmat, options, test_subject_projs) - symm_obj.run() - s_vec = SVectors(zmat, options, zmat.variable_dictionary_b) - s_vec.run(zmat.cartesians_b, True, second_order=options.second_order) - symm_obj.make_proj(s_vec) - assert np.allclose( - symm_obj.salc_proj, test_subject_projs[t], rtol=0, atol=1e-8 - ) - os.chdir("../../") diff --git a/concordantmodes/tests/nohup.out b/concordantmodes/tests/nohup.out deleted file mode 100644 index 7855dc5..0000000 --- a/concordantmodes/tests/nohup.out +++ /dev/null @@ -1,26 +0,0 @@ -============================= test session starts ============================== -platform linux -- Python 3.13.12, pytest-9.0.2, pluggy-1.6.0 -rootdir: /home/mel64643/github/concordantmodes -configfile: pyproject.toml -plugins: xdist-3.8.0, anyio-4.12.1 -collected 40 items - -test_directory_tree.py . [ 2%] -test_f_convert.py .. [ 7%] -test_f_read.py . [ 10%] -test_g_matrix.py . [ 12%] -test_gf_method.py . [ 15%] -test_reap.py . [ 17%] -test_s_vectors.py . [ 20%] -test_transf_disp.py . [ 22%] -test_zmat.py ............................... [100%] - -=============================== warnings summary =============================== -../gf_method.py:61 -concordantmodes/tests/test_gf_method.py::test_gf_method -concordantmodes/tests/test_transf_disp.py::test_transf_disp - /home/mel64643/github/concordantmodes/concordantmodes/gf_method.py:61: ComplexWarning: Casting complex values to real discards the imaginary part - self.freq = self.freq.astype(float) - --- Docs: https://docs.pytest.org/en/stable/how-to/capture-warnings.html -======================== 40 passed, 3 warnings in 8.15s ======================== diff --git a/concordantmodes/tests/suite_execute.py b/concordantmodes/tests/suite_execute.py index 9c38370..0e3c9c0 100644 --- a/concordantmodes/tests/suite_execute.py +++ b/concordantmodes/tests/suite_execute.py @@ -56,7 +56,8 @@ def run(self): self.ZMAT, "internal", False, - self.TED_obj, + self.s_vec.proj, + # self.TED_obj.proj, self.options, ) self.f_conv.run() @@ -64,11 +65,16 @@ def run(self): self.F = self.f_conv.F self.G = np.dot(self.TED_obj.proj.T, np.dot(self.g_mat.G, self.TED_obj.proj)) self.GF = GFMethod( - self.G, self.F, self.ZMAT, self.TED_obj, self.options, self.symm_obj.symtext + self.G, + self.F, + self.ZMAT, + self.TED_obj, + self.options, + symtext=self.symm_obj.symtext, ) self.GF.run() - self.algo = Algorithm(len(self.GF.L), None, self.options, None) + self.algo = Algorithm(len(self.GF.L), None, self.options) self.algo.run() self.disps = TransfDisp( @@ -76,7 +82,7 @@ def run(self): self.ZMAT, self.GF.L, True, - self.TED_obj, + self.TED_obj.proj, self.options, self.algo.indices, symm_obj=self.symm_obj, diff --git a/concordantmodes/tests/test_directory_tree.py b/concordantmodes/tests/test_directory_tree.py index 3c13d66..d69e25c 100644 --- a/concordantmodes/tests/test_directory_tree.py +++ b/concordantmodes/tests/test_directory_tree.py @@ -39,9 +39,15 @@ def test_make_input(): None, symm_obj, "template.dat", - None, + ".", + # None, ) + DT.init = False + DT.genbas = False + DT.ecp = False + DT.sub = False + with open("template.dat", "r") as file: data = file.readlines() with open("template_ref.dat", "r") as file: @@ -50,7 +56,8 @@ def test_make_input(): if os.path.exists(os.getcwd() + "/1"): shutil.rmtree(os.getcwd() + "/1") - data = DT.make_input(data, dispp, len(at), at, index, "input.dat", "1") + DT.make_input(data, dispp, len(at), "input.dat", "1") + # data = DT.make_input(data, dispp, len(at), at, index, "input.dat", "1") os.chdir("../..") - assert data == reference + assert DT.new_data == reference diff --git a/concordantmodes/tests/test_f_convert.py b/concordantmodes/tests/test_f_convert.py index 0fa40e9..424ce3f 100644 --- a/concordantmodes/tests/test_f_convert.py +++ b/concordantmodes/tests/test_f_convert.py @@ -28,7 +28,7 @@ def test_f_convert2int(): suite.ZMAT, "internal", False, - suite.TED_obj, + suite.TED_obj.proj, suite.options, ) FCint.run() @@ -49,7 +49,7 @@ def test_f_convert2cart(): suite.ZMAT, "internal", False, - suite.TED_obj, + suite.TED_obj.proj, suite.options, ) FCint.run() @@ -59,7 +59,7 @@ def test_f_convert2cart(): suite.ZMAT, "cartesian", False, - suite.TED_obj, + suite.TED_obj.proj, suite.options, ) FCcart.run() diff --git a/concordantmodes/tests/test_s_vectors.py b/concordantmodes/tests/test_s_vectors.py index fd56add..d0b1cec 100644 --- a/concordantmodes/tests/test_s_vectors.py +++ b/concordantmodes/tests/test_s_vectors.py @@ -12,6 +12,7 @@ from concordantmodes.zmat import Zmat suite = execute_suite("./ref_data/s_vec_test/", "Custom", s_vec_bool=True) +print(os.getcwd()) suite.run() # B-tensor ran on INTDER diff --git a/concordantmodes/tests/test_symmetry.py b/concordantmodes/tests/test_symmetry.py new file mode 100644 index 0000000..bde0c3a --- /dev/null +++ b/concordantmodes/tests/test_symmetry.py @@ -0,0 +1,172 @@ +import os +from concordantmodes.symmetry import Symmetry +from concordantmodes.options import Options +from concordantmodes.zmat import Zmat +from concordantmodes.s_vectors import SVectors + +# import molsym +import numpy as np + + +def test_symmetry(): + os.chdir("./ref_data/symmetry_test/") + test_subjects_del = ["c2v_water_zmat_del", "c3v_ammonia_zmat_del"] + test_subjects_custom = ["c2v_water_zmat_custom", "c3v_ammonia_zmat_custom"] + # delocalized first + test_subject_projs_del = [ + np.array( + [ + [-0.56086644, 0.43061449, -0.70710678], + [-0.56086644, 0.43061449, 0.70710678], + [0.60898085, 0.79318492, 0.0], + ] + ), + np.array( + [ + [ + -0.14852746, + -0.55791839, + -0.3291588, + -0.57011977, + -0.24149772, + 0.41828633, + ], + [ + -0.14852746, + -0.55791839, + -0.3291588, + 0.57011977, + -0.24149772, + -0.41828633, + ], + [-0.14852746, -0.55791839, 0.6583176, 0.0, 0.48299545, -0.0], + [0.23900037, -0.06362601, 0.42112908, 0.0, -0.5739944, 0.0], + [ + 0.23900037, + -0.06362601, + -0.21056454, + 0.36470848, + 0.2869972, + 0.49709374, + ], + [ + 0.23900037, + -0.06362601, + -0.21056454, + -0.36470848, + 0.2869972, + -0.49709374, + ], + [0.50413466, -0.13420931, -0.23650561, 0.0, 0.3223546, 0.0], + [ + -0.50413466, + 0.13420931, + -0.1182528, + 0.20481986, + 0.1611773, + 0.27916727, + ], + [ + 0.50413466, + -0.13420931, + 0.1182528, + 0.20481986, + -0.1611773, + 0.27916727, + ], + ] + ), + ] + + # start of customs + test_subject_projs_custom = [ + np.array( + [ + [-0.56086644, 0.43061449, -0.70710678], + [-0.56086644, 0.43061449, 0.70710678], + [0.60898085, 0.79318492, 0.0], + ] + ), + np.array( + [ + [ + -0.14852746, + -0.55791839, + -0.3291588, + 0.57011977, + -0.24149772, + 0.41828633, + ], + [ + -0.14852746, + -0.55791839, + -0.3291588, + -0.57011977, + -0.24149772, + -0.41828633, + ], + [-0.14852746, -0.55791839, 0.6583176, -0.0, 0.48299545, -0.0], + [0.23900037, -0.06362601, 0.42112908, 0.0, -0.5739944, 0.0], + [ + 0.23900037, + -0.06362601, + -0.21056454, + 0.36470848, + 0.2869972, + -0.49709374, + ], + [ + 0.23900037, + -0.06362601, + -0.21056454, + -0.36470848, + 0.2869972, + 0.49709374, + ], + [ + 0.50413466, + -0.13420931, + 0.1182528, + -0.20481986, + -0.1611773, + 0.27916727, + ], + [ + 0.50413466, + -0.13420931, + 0.1182528, + 0.20481986, + -0.1611773, + -0.27916727, + ], + [0.50413466, -0.13420931, -0.23650561, 0.0, 0.3223546, -0.0], + ] + ), + ] + # iterate over delocalized and custom coordinates + for c, coord in enumerate(["Delocalized", "Custom"]): + if coord == "Delocalized": + test_subjects = test_subjects_del + test_subject_projs = test_subject_projs_del + else: + test_subjects = test_subjects_custom + test_subject_projs = test_subject_projs_custom + + # iterate over zmats to test symmetry projection on + for t, test in enumerate(test_subjects): + options = Options() + options.cart_insert = 9 + options.molsym_symmetry = True + options.coords = coord + options.second_order = False + zmat = Zmat(options) + zmat.run(test) + symm_obj = Symmetry(zmat, options, test_subject_projs) + symm_obj.run() + s_vec = SVectors(zmat, options) + s_vec.run(zmat.cartesians_b, True, second_order=options.second_order) + symm_obj.make_proj(s_vec) + assert np.allclose( + symm_obj.salc_proj, test_subject_projs[t], rtol=0, atol=1e-8 + ) + os.chdir("../../") diff --git a/concordantmodes/tests/test_transf_disp.py b/concordantmodes/tests/test_transf_disp.py index a806989..7fdbbe2 100644 --- a/concordantmodes/tests/test_transf_disp.py +++ b/concordantmodes/tests/test_transf_disp.py @@ -31,7 +31,7 @@ def test_transf_disp(): suite.ZMAT, suite.GF.L, True, - suite.TED_obj, + suite.TED_obj.proj, suite.options, suite.algo.indices, symm_obj=suite.symm_obj, diff --git a/concordantmodes/transf_disp.py b/concordantmodes/transf_disp.py index 6eb7dad..19732ea 100644 --- a/concordantmodes/transf_disp.py +++ b/concordantmodes/transf_disp.py @@ -9,35 +9,131 @@ np.set_printoptions(precision=8, linewidth=240) -class TransfDisp(object): - np.set_printoptions(precision=8, linewidth=240) +class TransfDisp: """ - The purpose of this script is to perform an iterative transformation - between internal coordinate displacements and cartesian displacements, - using a first order, linear B-tensor transformation. - - It is necessary to construct an inverse B-tensor called 'A', where, - - A = uB^T(BuB^T)^-1. - - u can simply be the I_3N identity matrix, so long as frequency intensities - are not desired. - - If this is the case, then A simply becomes the following: - - A = B^T(BB^T)^-1. - - Where in fact, BB^T = G (non-mass-weighted) - + Generate Cartesian displacement geometries from internal-coordinate or + Cartesian-coordinate perturbations for finite-difference derivative + calculations. + + This class performs iterative transformations between normal-coordinate + displacements and Cartesian coordinates using Wilson's B-matrix formalism. + It constructs the generalized inverse transformation matrix (A-matrix), + generates displaced molecular geometries, and optionally includes + second-order corrections through numerical second derivatives of the + B-matrix. + + The transformation is based on: + + A = B⁺ P + + where B⁺ is the pseudoinverse of the B-matrix and P is the projection + matrix defining the nonredundant internal coordinate space. The resulting + A-matrix is subsequently transformed into the normal-coordinate basis to + produce Cartesian displacements corresponding to specified vibrational + normal mode displacements. + + The class supports: + + * Internal-coordinate displacements in normal-coordinate space. + * Cartesian-coordinate finite-difference displacements. + * First-derivative (gradient) displacement generation. + * Second-derivative (Hessian) displacement generation. + * Reduced or scaled displacement schemes. + * Optional second-order coordinate transformations using A₂ tensors. + * Tight iterative coordinate transformations with B-matrix updates. + * Symmetry-adapted displacement generation. + + Parameters + ---------- + s_vectors : SVectors or None + SVectors object containing the B-matrix and optional second-order + B-tensors. May be None for pure Cartesian displacements. + zmat : Zmat + Molecular coordinate object containing Cartesian coordinates, + internal-coordinate definitions, masses, and connectivity data. + eigs : ndarray + Matrix of normal-mode eigenvectors. + conv : bool + If True, enforce angular continuity corrections when converting + Cartesian coordinates to internal coordinates. + proj : ndarray + Projection matrix defining the nonredundant internal coordinate space. + options : Options + Runtime options controlling displacement sizes, convergence criteria, + scaling methods, coordinate systems, and numerical differentiation. + indices : iterable + List of coordinate index pairs used when generating second-derivative + displacements. + symm_obj : object, optional + Symmetry object used for symmetry-adapted coordinate generation. + coord_type : {"internal", "cartesian"}, optional + Coordinate system used for displacement generation. + Default is "internal". + deriv_level : int, optional + Derivative level to generate: + + * 0 : second derivatives (Hessian displacements) + * 1 : first derivatives (gradient displacements) + + cma_level : str, optional + Coordinate mode analysis level used when scaling displacements. + Default is "B". + + Attributes + ---------- + B : ndarray + Wilson B-matrix relating internal and Cartesian coordinates. + A : ndarray + First-order transformation matrix from normal coordinates to + Cartesian displacements. + eig_inv : ndarray + Normalized inverse eigenvector matrix. + n_coord : ndarray + Reference normal-coordinate values. + p_disp : ndarray + Positive displacement geometries. + m_disp : ndarray + Negative displacement geometries. + ref_carts : ndarray + Reference Cartesian geometry. + disp : float or ndarray + Displacement magnitudes used during finite differences. + proj : ndarray + Internal-coordinate projection matrix. + coord_type : str + Coordinate system used for displacement generation. + deriv_level : int + Derivative order being generated. + + Notes + ----- + For internal-coordinate displacements, Cartesian geometries are generated + iteratively until the transformed geometry reproduces the target normal + coordinate displacement within the specified convergence tolerance. + + When ``options.second_order`` is enabled, second-order coordinate + transformations are included through: + + Δx = A ΔQ + 1/2 A₂ : ΔQΔQ + + where A₂ is constructed from numerical second derivatives of the + Wilson B-matrix. + + The ordering of internal coordinates must remain consistent with the + ordering defined throughout the Concordant Modes package, since all + coordinate transformations and force-constant manipulations depend on + that convention. """ + np.set_printoptions(precision=8, linewidth=240) + def __init__( self, s_vectors, zmat, eigs, conv, - ted, + proj, options, indices, symm_obj=None, @@ -57,7 +153,7 @@ def __init__( self.ref_carts = np.array(self.ref_carts).astype(float) self.u = np.identity(3 * len(zmat.atom_list)) self.disp = self.options.disp - self.ted = ted + self.proj = proj self.eigs = eigs self.disp_cart = {} self.disp_cart["ref"] = self.ref_carts.copy() @@ -67,694 +163,396 @@ def __init__( self.coord_type = coord_type self.cma_level = cma_level - def run(self, fc=np.array([])): + def run(self, fc=None): np.set_printoptions(precision=8, linewidth=240) - # Invert the L-matrix and then normalize the rows. - proj_tol = 1.0e-3 - self.eig_inv = LA.inv(self.eigs) # (Normal modes (Q) x Proj internals (S) ) - for i in range(len(self.eig_inv)): - self.eig_inv[i] = self.eig_inv[i] / LA.norm(self.eig_inv[i]) - self.eig_inv[i][ - np.abs(self.eig_inv[i]) < np.max(np.abs(self.eig_inv[i])) * proj_tol - ] = 0 - - # Compute the A-matrix to convert from normal coordinates to cartesians - - u = np.array(self.zmat.masses) - for i in range(len(u)): - if self.zmat.atom_list[i] == "X": - u[i] = 0 - else: - u[i] = 1.0 / u[i] - u = np.repeat(u, 3) - u = np.diag(u) - # print(self.coord_type) - # raise RuntimeError - if self.coord_type != "cartesian": - # print(self.coord_type) - self.A = self.compute_A( - self.B.copy(), self.ted.proj, self.eig_inv, u, cma_level=self.cma_level - ) - - # Generate the normal coordinate values for the reference structure - self.n_coord = self.int_c(self.ref_carts, self.eig_inv, self.ted.proj) - - print("Normal Coordinate Values:") - for i in range(len(self.n_coord)): - print( - "Normal Coordinate #{:<4n}: {: 3.5f}".format(i + 1, self.n_coord[i]) - ) + fc = np.asarray(fc) if fc is not None else np.array([]) - # Next, we will have to specify our Normal mode internal coordinate - # displacement sizes + self._build_eig_inv() if self.coord_type == "internal": - self.Disp = self.disp - self.disp = [] - for i in range(len(self.n_coord)): - self.disp.append(self.Disp) - - # Reduced displacements generator - if self.options.reduced_disp and len(fc): - disp_size = self.options.reduced_disp_size - for i in range(len(self.disp)): - self.disp[i] = disp_size / np.abs(fc[i, i]) ** 0.25 - print("Reduced displacements") - print(self.disp) - buff_disp = np.zeros(len(self.disp)) - # Use this loop to print out the reduced disps in internal coords. - # Might want to make this an optional printout. - for i in range(len(self.disp)): - buff_disp[i] = self.disp[i] - simple_disp = np.dot(buff_disp, self.eig_inv) - simple_disp = np.dot(simple_disp, self.ted.proj.T) - print("Reduced displacement in simple internal coordinates:") - print(self.disp.copy()[i]) - print(simple_disp) - buff_disp *= 0 - - # This code loops through a list of indices that the force constants will be computed at and - # generates the displacements for the diagonal and (if specified) the off-diagonals. Where - # the displacement matrix D[i,j] = D[j,i]. - - disp = np.zeros(len(self.eigs.T)) - buff = disp.copy() - if not self.deriv_level: - p_disp = np.zeros((len(self.eigs), len(self.eigs)), dtype=object) - m_disp = np.zeros((len(self.eigs), len(self.eigs)), dtype=object) - disp = np.zeros((len(self.disp))) - if self.symm_obj.symtext is not None and self.options.exploit_pm_symm: - print( - "Molsym is being used. +/- displacements of non TSIR are equivalent." - ) - if self.options.only_TSIR: - indices = self.symm_obj.indices_by_irrep[0] - else: - print( - "Generating displacements for all irreps. Only + displacements for non TSIR" - ) - indices = ( - copy.deepcopy(self.symm_obj.indices_by_irrep) - .flatten() - .reshape((-1, 2)) - ) - else: - indices = self.indices - - Sum = 1 - h = 0 - for index in indices: - i, j = index[0], index[1] - disp[i] = self.disp[i] - disp[j] = self.disp[j] - p_disp[i, j] = self.coord_convert( - disp, - self.n_coord.copy(), - self.ref_carts.copy(), - self.A.copy(), - False, - self.zmat, - self.options, - ) - if h != 0: - pass - else: - m_disp[i, j] = self.coord_convert( - -disp, - self.n_coord.copy(), - self.ref_carts.copy(), - self.A.copy(), - False, - self.zmat, - self.options, - ) - disp = buff.copy() - if ( - self.symm_obj.symtext is not None - and self.options.exploit_pm_symm - and not self.options.only_TSIR - and Sum > len(self.symm_obj.indices_by_irrep[0]) - ): - h += 1 - Sum += 1 - elif self.deriv_level == 1: - p_disp = np.zeros(len(self.eigs), dtype=object) - m_disp = np.zeros(len(self.eigs), dtype=object) - for i in range(len(self.eigs)): - disp[i] = self.disp[i] - p_disp[i] = self.coord_convert( - disp, - self.n_coord.copy(), - self.ref_carts.copy(), - self.A.copy(), - False, - self.zmat, - self.options, - ) - m_disp[i] = self.coord_convert( - -disp, - self.n_coord.copy(), - self.ref_carts.copy(), - self.A.copy(), - False, - self.zmat, - self.options, - ) - disp = buff.copy() - else: - print( - "Only energy and gradient derivatives are supported. Check your deriv_level_b keyword." - ) - raise RuntimeError - self.p_disp = p_disp - self.m_disp = m_disp + self._run_internal(fc) elif self.coord_type == "cartesian": - self.disp_mag = self.options.disp - if self.deriv_level == 1: - p_disp = np.zeros(len(self.ref_carts.flatten()), dtype=object) - m_disp = np.zeros(len(self.ref_carts.flatten()), dtype=object) - for i in range(len(self.ref_carts.flatten())): - p_disp[i] = self.ref_carts.flatten().copy() - m_disp[i] = self.ref_carts.flatten().copy() - p_disp[i][i] = p_disp[i][i] + self.disp - m_disp[i][i] = m_disp[i][i] - self.disp - p_disp[i] = np.reshape(p_disp[i], (-1, 3)) - m_disp[i] = np.reshape(m_disp[i], (-1, 3)) - - self.p_disp = p_disp - self.m_disp = m_disp - elif not self.deriv_level: - if self.options.molsym_symmetry: - print("The ref geom") - print(self.ref_carts) - print(self.symm_obj.symtext.mol.coords) - p_disp = np.zeros( - (len(self.ref_carts.flatten()), len(self.ref_carts.flatten())), - dtype=object, - ) - m_disp = np.zeros( - (len(self.ref_carts.flatten()), len(self.ref_carts.flatten())), - dtype=object, - ) - print(f"p disp {p_disp}") - n_irrep = len(self.symm_obj.symtext.irreps) - salc_indices_pi = self.symm_obj.CDsalcs.salcs_by_irrep - print("salc indices pi") - print(salc_indices_pi) - new_indices = [] - for i in range(len(self.symm_obj.CDsalcs)): - for j in range(len(self.symm_obj.CDsalcs)): - print(f"i j pair {i,j}") - pair = [i, j] - disp_geom = self.displace_geom(pair) - print("the disp geom") - print(disp_geom) - if i == j: - p_disp[i, i] += disp_geom - m_disp[i, i] -= disp_geom - else: - p_disp[i, j] += disp_geom - m_disp[i, j] -= disp_geom - # print(stop) - - # for i in range(len(self.symm_obj.CDsalcs)): - # for j in range(len(self.symm_obj.CDsalcs)): - # if i <= j: - # salc_i = self.symm_obj.CDsalcs[i] - # salc_j = self.symm_obj.CDsalcs[j] - # print(f"the salc {i, j}") - # print(f"The disp {self.disp}") - # new_indices.append([i,j]) - # #self.append_geoms(salc_i, salc_j) - # disp_geom = np.copy(self.symm_obj.symtext.mol.coords) - # for atom_idx in range(self.symm_obj.symtext.mol.natoms): - # for xyz_idx in range(3): - # disp_geom[atom_idx, xyz_idx] += self.disp * salc_i.coeffs[atom_idx * 3 + xyz_idx] #/ np.sqrt(self.symm_obj.symtext.mol.masses[atom_idx]) - # disp_geom[atom_idx, xyz_idx] += self.disp * salc_j.coeffs[atom_idx * 3 + xyz_idx] #/ np.sqrt(self.symm_obj.symtext.mol.masses[atom_idx]) - # print(f"The disp_geom {disp_geom}") - # if i == j: - # #p_disp[i,i][i] += self.disp - # #m_disp[i,i][i] -= self.disp - # p_disp[i,i] += disp_geom - # m_disp[i,i] -= disp_geom - # else: - # #p_disp[i,j][i] += self.disp - # #p_disp[i,j][j] += self.disp - # #m_disp[i,j][i] -= self.disp - # #m_disp[i,j][j] -= self.disp - # p_disp[i,j] += disp_geom - # m_disp[i,j] -= disp_geom - self.p_disp = p_disp - self.m_disp = m_disp - else: - p_disp = np.zeros( - (len(self.ref_carts.flatten()), len(self.ref_carts.flatten())), - dtype=object, - ) - m_disp = np.zeros( - (len(self.ref_carts.flatten()), len(self.ref_carts.flatten())), - dtype=object, - ) - # print(self.indices) - # raise RuntimeError - for index in self.indices: - i, j = index[0], index[1] - if i == j: - p_disp[i, i] = self.ref_carts.flatten().copy() - m_disp[i, i] = self.ref_carts.flatten().copy() - p_disp[i, i][i] += self.disp - m_disp[i, i][i] -= self.disp - p_disp[i, i] = np.reshape(p_disp[i, i], (-1, 3)) - m_disp[i, i] = np.reshape(m_disp[i, i], (-1, 3)) - else: - p_disp[i, j] = self.ref_carts.flatten().copy() - m_disp[i, j] = self.ref_carts.flatten().copy() - p_disp[i, j][i] += self.disp - p_disp[i, j][j] += self.disp - m_disp[i, j][i] -= self.disp - m_disp[i, j][j] -= self.disp - p_disp[i, j] = np.reshape(p_disp[i, j], (-1, 3)) - m_disp[i, j] = np.reshape(m_disp[i, j], (-1, 3)) - - self.p_disp = p_disp - self.m_disp = m_disp + self._run_cartesian() else: - print( - "Please input a displacement coordinate type of either cartesian or internal." + raise RuntimeError("coord_type must be either 'cartesian' or 'internal'.") + + def _build_eig_inv(self, proj_tol=1.0e-3): + """ + Invert and normalize the eigenvector matrix. + """ + + self.eig_inv = LA.inv(self.eigs) + + for i, row in enumerate(self.eig_inv): + + row /= LA.norm(row) + + thresh = np.max(np.abs(row)) * proj_tol + row[np.abs(row) < thresh] = 0.0 + + self.eig_inv[i] = row + + def _build_mass_matrix(self): + """ + Construct inverse mass-weight matrix. + """ + + masses = np.asarray(self.zmat.masses, dtype=float) + + inv_masses = np.where( + np.asarray(self.zmat.atom_list) == "X", + 0.0, + 1.0 / masses, + ) + + return np.diag(np.repeat(inv_masses, 3)) + + # Internal displacement functions here: + def _run_internal(self, fc): + + u = self._build_mass_matrix() + + self.A = self.compute_A( + self.B.copy(), + self.proj, + self.eig_inv, + u, + cma_level=self.cma_level, + ) + + self.n_coord = self.int_c( + self.ref_carts, + self.eig_inv, + self.proj, + ) + + self._print_normal_coordinates() + + self._build_displacements(fc) + + A2 = self._build_second_order_A2() + + if self.deriv_level == 1: + self._generate_internal_first_deriv(A2) + + elif self.deriv_level == 0: + self._generate_internal_second_deriv(A2) + + else: + raise RuntimeError("Only energy and gradient derivatives are supported.") + + def _print_normal_coordinates(self): + + print("Normal Coordinate Values:") + + for i, value in enumerate(self.n_coord, start=1): + print(f"Normal Coordinate #{i:<4}: {value: 3.5f}") + + def _build_displacements(self, fc): + + self.Disp = self.disp + self.disp = np.full(len(self.n_coord), self.Disp) + + # + # Reduced displacements + # + if self.options.reduced_disp and len(fc): + + scale = self.options.reduced_disp_size + + self.disp = np.array( + [scale / abs(fc[i, i]) ** 0.25 for i in range(len(self.disp))] + ) + + print("Reduced displacements") + print(self.disp) + + # + # CMA scaling + # + elif self.options.scaled_disp and self.cma_level == "B": + + for i in range(len(self.disp)): + + self.disp[i] /= np.max(self.proj.T[i]) + self.disp[i] /= np.max(self.eig_inv[i]) + self.disp[i] *= LA.norm(self.eig_inv[i]) + + def _build_second_order_A2(self): + + if not self.options.second_order: + return None + + L = inv(self.eig_inv) + + B2 = np.einsum( + "ir,ipq->rpq", + self.proj, + self.s_vectors.B2, + ) + + B2 = np.einsum( + "ri,ipq->rpq", + L, + B2, + ) + + return self.compute_A2(B2, self.A) + + def _generate_internal_first_deriv(self, A2): + + n = len(self.eigs) + + p_disp = np.empty(n, dtype=object) + m_disp = np.empty(n, dtype=object) + + for i in range(n): + + disp = np.zeros(n) + disp[i] = self.disp[i] + + p_disp[i] = self.coord_convert( + disp, + self.n_coord.copy(), + self.ref_carts.copy(), + self.A.copy(), + False, + self.zmat, + self.options, + A2=A2, + ) + + m_disp[i] = self.coord_convert( + -disp, + self.n_coord.copy(), + self.ref_carts.copy(), + self.A.copy(), + False, + self.zmat, + self.options, + A2=A2, ) - raise RuntimeError - - # function useful for CMA geometry optimizer - def eigs_inv(self): - proj_tol = 1.0e-3 - self.eig_inv = LA.inv(self.eigs) # (Normal modes (Q) x Sym internals (S) ) - for i in range(len(self.eig_inv)): - self.eig_inv[i] = self.eig_inv[i] / LA.norm(self.eig_inv[i]) - self.eig_inv[i][ - np.abs(self.eig_inv[i]) < np.max(np.abs(self.eig_inv[i])) * proj_tol - ] = 0 - - def displace_geom(self, pair): - print("inside displace geom") - print(f"The pair {pair}") - disp_geom = np.copy(self.symm_obj.symtext.mol.coords) - for s_index in pair: - print(f"The s index {s_index}") - for atom_idx in range(self.symm_obj.symtext.mol.natoms): - for xyz_idx in range(3): - disp_geom[atom_idx, xyz_idx] += ( - self.disp - * self.symm_obj.CDsalcs.salcs[s_index].coeffs[ - atom_idx * 3 + xyz_idx - ] - ) # / np.sqrt(self.symm_obj.symtext.mol.masses[atom_idx]) - return disp_geom + + self.p_disp = p_disp + self.m_disp = m_disp + + def _generate_internal_second_deriv(self, A2): + + n = len(self.eigs) + + p_disp = np.empty((n, n), dtype=object) + m_disp = np.empty((n, n), dtype=object) + + for i, j in self.indices: + + disp = np.zeros(n) + + disp[i] = self.disp[i] + disp[j] = self.disp[j] + + p_disp[i, j] = self.coord_convert( + disp, + self.n_coord.copy(), + self.ref_carts.copy(), + self.A.copy(), + self.options.tight_disp, + self.zmat, + self.options, + A2=A2, + ) + + m_disp[i, j] = self.coord_convert( + -disp, + self.n_coord.copy(), + self.ref_carts.copy(), + self.A.copy(), + self.options.tight_disp, + self.zmat, + self.options, + A2=A2, + ) + + self.p_disp = p_disp + self.m_disp = m_disp + + # Now we have the cartesian functions. + def _run_cartesian(self): + + self.disp_mag = self.options.disp + if self.deriv_level == 1: + self._generate_cartesian_first_deriv() + + elif self.deriv_level == 0: + + if not self.options.molsym_symmetry: + self._generate_cartesian_second_deriv() + + def _generate_cartesian_first_deriv(self): + + ref = self.ref_carts.flatten() + n = len(ref) + + p_disp = np.empty(n, dtype=object) + m_disp = np.empty(n, dtype=object) + + for i in range(n): + + plus = ref.copy() + minus = ref.copy() + + plus[i] += self.disp + minus[i] -= self.disp + + p_disp[i] = plus.reshape(-1, 3) + m_disp[i] = minus.reshape(-1, 3) + + self.p_disp = p_disp + self.m_disp = m_disp + + def _generate_cartesian_second_deriv(self): + + ref = self.ref_carts.flatten() + n = len(ref) + + p_disp = np.empty((n, n), dtype=object) + m_disp = np.empty((n, n), dtype=object) + + for i, j in self.indices: + + plus = ref.copy() + minus = ref.copy() + + plus[i] += self.disp + plus[j] += self.disp + + minus[i] -= self.disp + minus[j] -= self.disp + + p_disp[i, j] = plus.reshape(-1, 3) + m_disp[i, j] = minus.reshape(-1, 3) + + self.p_disp = p_disp + self.m_disp = m_disp def int_c(self, carts, eig_inv, proj): - # This is a function that computes all currently implemented and - # specified internal coordinates from the desired cartesian - # coordinates. The internal coordinate values will be appended - # together in the order of the for-loops below. All other methods - # in this package which use the int_c generated variables MUST - # abide by this ordering. - int_coord = np.array([]) - for i in range(len(self.zmat.bond_indices)): - for j in self.zmat.bond_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.bond_indices[i] - bond_carts = [] - for k in indies[0]: - bond_carts = np.append( - bond_carts, np.array(carts[int(k) - 1]).astype(float) - ) - bond_carts = np.reshape(bond_carts, (-1, 3)) - x1 = self.calc_Centroid(bond_carts) - - bond_carts = [] - for k in indies[1]: - bond_carts = np.append( - bond_carts, np.array(carts[int(k) - 1]).astype(float) - ) - bond_carts = np.reshape(bond_carts, (-1, 3)) - x2 = self.calc_Centroid(bond_carts) - else: - x1 = np.array(carts[int(self.zmat.bond_indices[i][0]) - 1]).astype( - float - ) - x2 = np.array(carts[int(self.zmat.bond_indices[i][1]) - 1]).astype( - float - ) - - int_coord = np.append(int_coord, self.calc_bond(x1, x2)) - - for i in range(len(self.zmat.rcom_indices)): - mol1 = self.zmat.rcom_indices[i][0] - mol2 = self.zmat.rcom_indices[i][1] - rc = self.calc_Rcom(mol1, mol2, carts) - int_coord = np.append(int_coord, rc) - - for i in range(len(self.zmat.angle_indices)): - for j in self.zmat.angle_indices[i]: - Len = len(np.array(j).shape) - indies = self.zmat.angle_indices[i] - if Len: - angle_carts = [] - for k in indies[0]: - angle_carts = np.append( - angle_carts, np.array(carts[int(k) - 1]).astype(float) - ) - angle_carts = np.reshape(angle_carts, (-1, 3)) - x1 = self.calc_Centroid(angle_carts) - - angle_carts = [] - for k in indies[1]: - angle_carts = np.append( - angle_carts, np.array(carts[int(k) - 1]).astype(float) - ) - angle_carts = np.reshape(angle_carts, (-1, 3)) - x2 = self.calc_Centroid(angle_carts) - - angle_carts = [] - for k in indies[2]: - angle_carts = np.append( - angle_carts, np.array(carts[int(k) - 1]).astype(float) - ) - angle_carts = np.reshape(angle_carts, (-1, 3)) - x3 = self.calc_Centroid(angle_carts) - else: - x1 = np.array(carts[int(indies[0]) - 1]).astype(float) - x2 = np.array(carts[int(indies[1]) - 1]).astype(float) - x3 = np.array(carts[int(indies[2]) - 1]).astype(float) - - a = self.calc_angle(x1, x2, x3) + """ + Compute all internal coordinates from Cartesian coordinates. + + The ordering of generated internal coordinates is significant and must + remain consistent with the rest of the package. + """ + carts = np.asarray(carts, dtype=float) + int_coord = [] + + def get_point(indices): + """ + Return either: + - a Cartesian coordinate for a single atom index + - a centroid for a group of atom indices + """ + if np.ndim(indices) == 0: + return carts[int(indices) - 1] + + coords = np.asarray( + [carts[int(i) - 1] for i in indices], + dtype=float, + ) + return self.calc_Centroid(coords) + + def get_points(index_set): + return [get_point(indices) for indices in index_set] + + # + # Bonds + # + for inds in self.zmat.bond_indices: + x1, x2 = get_points(inds) + int_coord.append(self.calc_bond(x1, x2)) + + # + # Center-of-mass distances + # + for mol1, mol2 in self.zmat.rcom_indices: + int_coord.append(self.calc_Rcom(mol1, mol2, carts)) + + # + # Angles + # + for i, inds in enumerate(self.zmat.angle_indices): + x1, x2, x3 = get_points(inds) + + angle = self.calc_angle(x1, x2, x3) + if self.conv: - condition_1 = ( - float(self.zmat.variable_dictionary_a[self.zmat.angle_variables[i]]) - > 180.0 - ) - condition_2 = ( - float(self.zmat.variable_dictionary_a[self.zmat.angle_variables[i]]) - < 0.0 + ref = float( + self.zmat.variable_dictionary_a[self.zmat.angle_variables[i]] ) - if condition_1 or condition_2: - a = 2 * np.pi - a - int_coord = np.append(int_coord, a) - - for i in range(len(self.zmat.torsion_indices)): - for j in self.zmat.torsion_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.torsion_indices[i] - torsion_carts = [] - for k in indies[0]: - torsion_carts = np.append( - torsion_carts, np.array(carts[int(k) - 1]).astype(float) - ) - torsion_carts = np.reshape(torsion_carts, (-1, 3)) - x1 = self.calc_Centroid(torsion_carts) - - torsion_carts = [] - for k in indies[1]: - torsion_carts = np.append( - torsion_carts, np.array(carts[int(k) - 1]).astype(float) - ) - torsion_carts = np.reshape(torsion_carts, (-1, 3)) - x2 = self.calc_Centroid(torsion_carts) - - torsion_carts = [] - for k in indies[2]: - torsion_carts = np.append( - torsion_carts, np.array(carts[int(k) - 1]).astype(float) - ) - torsion_carts = np.reshape(torsion_carts, (-1, 3)) - x3 = self.calc_Centroid(torsion_carts) - - torsion_carts = [] - for k in indies[3]: - torsion_carts = np.append( - torsion_carts, np.array(carts[int(k) - 1]).astype(float) - ) - torsion_carts = np.reshape(torsion_carts, (-1, 3)) - x4 = self.calc_Centroid(torsion_carts) - else: - x1 = np.array( - carts[int(self.zmat.torsion_indices[i][0]) - 1] - ).astype(float) - x2 = np.array( - carts[int(self.zmat.torsion_indices[i][1]) - 1] - ).astype(float) - x3 = np.array( - carts[int(self.zmat.torsion_indices[i][2]) - 1] - ).astype(float) - x4 = np.array( - carts[int(self.zmat.torsion_indices[i][3]) - 1] - ).astype(float) - - t = self.calc_tors(x1, x2, x3, x4) + + if ref > 180.0 or ref < 0.0: + angle = 2 * np.pi - angle + + int_coord.append(angle) + + # + # Torsions + # + for i, inds in enumerate(self.zmat.torsion_indices): + x1, x2, x3, x4 = get_points(inds) + + tors = self.calc_tors(x1, x2, x3, x4) + if self.conv: - condition_1 = float( - self.zmat.variable_dictionary_a[self.zmat.torsion_variables[i]] - ) > 135.0 and (t * 180.0 / np.pi < -135.0) - condition_2 = float( + ref = float( self.zmat.variable_dictionary_a[self.zmat.torsion_variables[i]] - ) < -135.0 and (t * 180.0 / np.pi > 135.0) - if condition_1: - t += 2 * np.pi - if condition_2: - t -= 2 * np.pi - int_coord = np.append(int_coord, t) - - for i in range(len(self.zmat.oop_indices)): - for j in self.zmat.oop_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.oop_indices[i] - oop_carts = [] - for k in indies[0]: - oop_carts = np.append( - oop_carts, np.array(carts[int(k) - 1]).astype(float) - ) - oop_carts = np.reshape(oop_carts, (-1, 3)) - x1 = self.calc_Centroid(oop_carts) - - oop_carts = [] - for k in indies[1]: - oop_carts = np.append( - oop_carts, np.array(carts[int(k) - 1]).astype(float) - ) - oop_carts = np.reshape(oop_carts, (-1, 3)) - x2 = self.calc_Centroid(oop_carts) - - oop_carts = [] - for k in indies[2]: - oop_carts = np.append( - oop_carts, np.array(carts[int(k) - 1]).astype(float) - ) - oop_carts = np.reshape(oop_carts, (-1, 3)) - x3 = self.calc_Centroid(oop_carts) - - oop_carts = [] - for k in indies[3]: - oop_carts = np.append( - oop_carts, np.array(carts[int(k) - 1]).astype(float) - ) - oop_carts = np.reshape(oop_carts, (-1, 3)) - x4 = self.calc_Centroid(oop_carts) - else: - x1 = np.array(carts[int(self.zmat.oop_indices[i][0]) - 1]).astype( - float - ) - x2 = np.array(carts[int(self.zmat.oop_indices[i][1]) - 1]).astype( - float - ) - x3 = np.array(carts[int(self.zmat.oop_indices[i][2]) - 1]).astype( - float - ) - x4 = np.array(carts[int(self.zmat.oop_indices[i][3]) - 1]).astype( - float - ) - - o = self.calc_OOP(x1, x2, x3, x4) - if self.conv: - condition_1 = ( - float(self.zmat.variable_dictionary_a[self.zmat.oop_variables[i]]) - > 180.0 - ) - condition_2 = ( - float(self.zmat.variable_dictionary_a[self.zmat.oop_variables[i]]) - < -180.0 ) - if condition_1: - o = o - 2 * np.pi - if condition_2: - o = o + 2 * np.pi - int_coord = np.append(int_coord, o) - - # These angles will have to be tempered at some point in the same manner as above. - for i in range(len(self.zmat.lin_indices)): - for j in self.zmat.lin_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.lin_indices[i] - lin_carts = [] - for k in indies[0]: - lin_carts = np.append( - lin_carts, np.array(carts[int(k) - 1]).astype(float) - ) - lin_carts = np.reshape(lin_carts, (-1, 3)) - x1 = self.calc_Centroid(lin_carts) - - lin_carts = [] - for k in indies[1]: - lin_carts = np.append( - lin_carts, np.array(carts[int(k) - 1]).astype(float) - ) - lin_carts = np.reshape(lin_carts, (-1, 3)) - x2 = self.calc_Centroid(lin_carts) - - lin_carts = [] - for k in indies[2]: - lin_carts = np.append( - lin_carts, np.array(carts[int(k) - 1]).astype(float) - ) - lin_carts = np.reshape(lin_carts, (-1, 3)) - x3 = self.calc_Centroid(lin_carts) - - lin_carts = [] - for k in indies[3]: - lin_carts = np.append( - lin_carts, np.array(carts[int(k) - 1]).astype(float) - ) - lin_carts = np.reshape(lin_carts, (-1, 3)) - x4 = self.calc_Centroid(lin_carts) - else: - x1 = np.array(carts[int(self.zmat.lin_indices[i][0]) - 1]).astype( - float - ) - x2 = np.array(carts[int(self.zmat.lin_indices[i][1]) - 1]).astype( - float - ) - x3 = np.array(carts[int(self.zmat.lin_indices[i][2]) - 1]).astype( - float - ) - x4 = np.array(carts[int(self.zmat.lin_indices[i][3]) - 1]).astype( - float - ) - l = self.calc_Lin(x1, x2, x3, x4) - int_coord = np.append(int_coord, l) - - for i in range(len(self.zmat.linx_indices)): - for j in self.zmat.linx_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.linx_indices[i] - linx_carts = [] - for k in indies[0]: - linx_carts = np.append( - linx_carts, np.array(carts[int(k) - 1]).astype(float) - ) - linx_carts = np.reshape(linx_carts, (-1, 3)) - x1 = self.calc_Centroid(linx_carts) - - linx_carts = [] - for k in indies[1]: - linx_carts = np.append( - linx_carts, np.array(carts[int(k) - 1]).astype(float) - ) - linx_carts = np.reshape(linx_carts, (-1, 3)) - x2 = self.calc_Centroid(linx_carts) - - linx_carts = [] - for k in indies[2]: - linx_carts = np.append( - linx_carts, np.array(carts[int(k) - 1]).astype(float) - ) - linx_carts = np.reshape(linx_carts, (-1, 3)) - x3 = self.calc_Centroid(linx_carts) - - linx_carts = [] - for k in indies[3]: - linx_carts = np.append( - linx_carts, np.array(carts[int(k) - 1]).astype(float) - ) - linx_carts = np.reshape(linx_carts, (-1, 3)) - x4 = self.calc_Centroid(linx_carts) - else: - x1 = np.array(carts[int(self.zmat.linx_indices[i][0]) - 1]).astype( - float - ) - x2 = np.array(carts[int(self.zmat.linx_indices[i][1]) - 1]).astype( - float - ) - x3 = np.array(carts[int(self.zmat.linx_indices[i][2]) - 1]).astype( - float - ) - x4 = np.array(carts[int(self.zmat.linx_indices[i][3]) - 1]).astype( - float - ) - lx = self.calc_Linx(x1, x2, x3, x4) - int_coord = np.append(int_coord, lx) - - for i in range(len(self.zmat.liny_indices)): - for j in self.zmat.liny_indices[i]: - Len = len(np.array(j).shape) - if Len: - indies = self.zmat.liny_indices[i] - liny_carts = [] - for k in indies[0]: - liny_carts = np.append( - liny_carts, np.array(carts[int(k) - 1]).astype(float) - ) - liny_carts = np.reshape(liny_carts, (-1, 3)) - x1 = self.calc_Centroid(liny_carts) - - liny_carts = [] - for k in indies[1]: - liny_carts = np.append( - liny_carts, np.array(carts[int(k) - 1]).astype(float) - ) - liny_carts = np.reshape(liny_carts, (-1, 3)) - x2 = self.calc_Centroid(liny_carts) - - liny_carts = [] - for k in indies[2]: - liny_carts = np.append( - liny_carts, np.array(carts[int(k) - 1]).astype(float) - ) - liny_carts = np.reshape(liny_carts, (-1, 3)) - x3 = self.calc_Centroid(liny_carts) - - liny_carts = [] - for k in indies[3]: - liny_carts = np.append( - liny_carts, np.array(carts[int(k) - 1]).astype(float) - ) - liny_carts = np.reshape(liny_carts, (-1, 3)) - x4 = self.calc_Centroid(liny_carts) - else: - x1 = np.array(carts[int(self.zmat.liny_indices[i][0]) - 1]).astype( - float - ) - x2 = np.array(carts[int(self.zmat.liny_indices[i][1]) - 1]).astype( - float - ) - x3 = np.array(carts[int(self.zmat.liny_indices[i][2]) - 1]).astype( - float - ) - x4 = np.array(carts[int(self.zmat.liny_indices[i][3]) - 1]).astype( - float - ) - ly = self.calc_Liny(x1, x2, x3, x4) - int_coord = np.append(int_coord, ly) - - int_coord = np.dot(proj.T, int_coord) - int_coord = np.dot(eig_inv, int_coord) - - return int_coord + + tors_deg = tors * 180.0 / np.pi + + if ref > 135.0 and tors_deg < -135.0: + tors += 2 * np.pi + elif ref < -135.0 and tors_deg > 135.0: + tors -= 2 * np.pi + + int_coord.append(tors) + + # + # Out-of-plane bends + # + for i, inds in enumerate(self.zmat.oop_indices): + x1, x2, x3, x4 = get_points(inds) + + oop = self.calc_OOP(x1, x2, x3, x4) + + if self.conv: + ref = float(self.zmat.variable_dictionary_a[self.zmat.oop_variables[i]]) + + if ref > 180.0: + oop -= 2 * np.pi + elif ref < -180.0: + oop += 2 * np.pi + + int_coord.append(oop) + + # + # Linear bends + # + linear_coordinate_sets = [ + (self.zmat.lin_indices, self.calc_Lin), + (self.zmat.linx_indices, self.calc_Linx), + (self.zmat.liny_indices, self.calc_Liny), + ] + + for index_sets, func in linear_coordinate_sets: + for inds in index_sets: + x1, x2, x3, x4 = get_points(inds) + int_coord.append(func(x1, x2, x3, x4)) + + int_coord = np.asarray(int_coord) + + return eig_inv @ (proj.T @ int_coord) def calc_bond(self, x1, x2): r = LA.norm(x1 - x2) @@ -868,12 +666,11 @@ def coord_convert( n_disp, n_coord, ref_carts, - # max_iter, - # tolerance, A, tight_disp, zmat, options, + A2=np.array([]), ): disp = n_disp.copy() new_n = n_coord + n_disp @@ -882,22 +679,18 @@ def coord_convert( for i in range(options.transf_max_it): cart_disp = np.dot(n_disp, A) np.set_printoptions(precision=6, linewidth=240) - # if self.options.second_order: - # A2 = self.compute_A2(self.s_vectors.B2, A) - # A2 = np.einsum("ipq,pr,qs->irs", A2, self.ted.proj, self.ted.proj) - # L = inv(self.eig_inv) - # A2 = np.einsum("irs,rp,sq->ipq",A2, L, L) - # cart_disp += 0.5*np.einsum("ipq,p,q->i", A2, n_disp,n_disp) + if self.options.second_order: + cart_disp += 0.5 * np.einsum("ipq,p,q->i", A2, n_disp, n_disp) cart_disp_shaped = np.reshape(cart_disp, (-1, 3)) new_carts += cart_disp_shaped - coord_check = self.int_c(new_carts, self.eig_inv, self.ted.proj) + coord_check = self.int_c(new_carts, self.eig_inv, self.proj) n_disp = new_n - coord_check if tight_disp: - sVec = s_vec(zmat, options, zmat.variable_dictionary_a) + sVec = s_vec(zmat, options) sVec.run(new_carts, False) A = self.compute_A( - sVec.B, self.ted.proj, self.eig_inv, self.zmat.mass_weight + sVec.B, self.proj, self.eig_inv, self.zmat.mass_weight ) if LA.norm(n_disp) < tolerance: break @@ -940,9 +733,8 @@ def compute_A(self, B, proj, eig_inv, u, cma_level="B"): return A def compute_A2(self, B2, A): + C2 = np.einsum("rij,pi,qj->rpq", B2, A, A) - C2 = np.einsum("rij,ip,jq->rpq", B2, A, A) - - A2 = np.einsum("rpq,ir->ipq", C2, A) + A2 = np.einsum("rpq,ri->ipq", C2, A) return A2 diff --git a/concordantmodes/vulcan_template.py b/concordantmodes/vulcan_template.py index 665a6ae..5052b7e 100644 --- a/concordantmodes/vulcan_template.py +++ b/concordantmodes/vulcan_template.py @@ -3,7 +3,7 @@ from numpy import linalg as LA -class VulcanTemplate(object): +class VulcanTemplate: """ This file just stores the vulcan qsub script """ diff --git a/concordantmodes/zmat.py b/concordantmodes/zmat.py index 39e4cc0..58aa0e9 100644 --- a/concordantmodes/zmat.py +++ b/concordantmodes/zmat.py @@ -9,7 +9,146 @@ from concordantmodes.transf_disp import TransfDisp -class Zmat(object): +class Zmat: + """ + Parse, construct, and analyze molecular internal coordinate systems. + + The ``Zmat`` class reads molecular geometry information from a ZMAT-style + input file, constructs internal coordinate definitions, evaluates their + numerical values from Cartesian coordinates, and provides mappings between + coordinate labels, atom indices, and coordinate values. + + The class supports three coordinate generation modes: + + * ``ZMAT`` – Traditional Z-matrix coordinates generated directly from + the atom connectivity specified in the Z-matrix section. + * ``DELOCALIZED`` – Redundant internal coordinates generated + automatically from molecular connectivity determined either from + user-supplied bonds or covalent radii analysis. + * ``CUSTOM`` – User-defined internal coordinates, including support + for centroid-based coordinates and center-of-mass distances. + + Supported internal coordinate types include: + + * Bond stretches + * Bond angles + * Dihedral/torsional angles + * Out-of-plane bends + * Linear bending coordinates + * Linear bend x/y components + * Center-of-mass separation coordinates + * Centroid-based coordinates + + The class reads one or two Cartesian geometries from the input file. + When two geometries are supplied, the first is treated as the initial + geometry and the second as the final geometry. Internal coordinate + values are evaluated for both geometries and stored for comparison. + + Attributes + ---------- + options : object + User options controlling coordinate generation, units, + topology analysis, and parsing behavior. + + atom_list : list[str] + Atomic labels read from the Cartesian geometry section. + + cartesians_b : ndarray, shape (N, 3) + Initial Cartesian coordinates in atomic units. + + cartesians_a : ndarray, shape (N, 3) + Final Cartesian coordinates in atomic units. + + masses : list[float] + Atomic masses converted to electron-mass units. + + mass_weight : ndarray, shape (3N, 3N) + Diagonal mass-weighting matrix. + + bond_indices : ndarray + Bond coordinate atom index definitions. + + angle_indices : ndarray + Angle coordinate atom index definitions. + + torsion_indices : ndarray + Torsional coordinate atom index definitions. + + oop_indices : ndarray + Out-of-plane bend coordinate definitions. + + lin_indices : ndarray + Linear bend coordinate definitions. + + linx_indices : ndarray + Linear bend x-component definitions. + + liny_indices : ndarray + Linear bend y-component definitions. + + rcom_indices : ndarray + Center-of-mass distance coordinate definitions. + + variables : ndarray + Ordered list of internal coordinate labels. + + variable_dictionary_b : dict + Internal coordinate values for the initial geometry. + + variable_dictionary_a : dict + Internal coordinate values for the final geometry. + + index_dictionary : dict + Mapping between coordinate labels and their corresponding + atom index definitions. + + Methods + ------- + run(zmat_name='zmat') + Execute the complete ZMAT processing workflow. + + zmat_read(zmat_name) + Read Cartesian coordinates and Z-matrix information from + an input file. + + zmat_process(zmat_output) + Construct internal coordinate definitions from parsed input. + + zmat_calc() + Compute internal coordinate values from Cartesian geometries. + + zmat_compile() + Build lookup dictionaries for coordinate labels and indices. + + zmat_print() + Print initial, final, and displacement values for all + internal coordinates. + + red_mass(indices) + Compute the reciprocal reduced mass contribution for a + coordinate definition. + + np_contains(array1, array2, tor=False) + Determine whether a coordinate definition already exists + within a coordinate list. + + Notes + ----- + Cartesian coordinates may be supplied in either Angstroms or + Bohr depending on the value of ``options.cart_coords``. All + internal calculations are performed in atomic units. + + For ``DELOCALIZED`` coordinates, molecular connectivity can be + determined automatically using covalent radii from + ``qcelemental.covalent_radii.CovalentRadii``. From this + connectivity graph, redundant bonds, angles, torsions, and + out-of-plane bends are generated automatically. + + When topology analysis is enabled, the class can additionally + generate and analyze molecular graph walks and ring cycles of + increasing length. + """ + def __init__(self, options): self.amu_elMass = 5.48579909065 * (10 ** (-4)) self.disp_tol = 1.0e-14 @@ -206,11 +345,6 @@ def zmat_process(self, zmat_output): # Second atom of the ZMAT, will have one bond term if re.search(self.second_atom_regex, zmat_output[i]): List = re.findall(self.second_atom_regex, zmat_output[i])[0] - print("Zebra2") - print(List) - # self.bond_indices = np.append( - # self.bond_indices, np.array([str(i - first_index + 1), List[1]]) - # ) self.bond_indices = np.append(self.bond_indices, 0) self.bond_indices[-1] = np.array( [str(i - first_index + 1), List[1]], dtype=object @@ -224,9 +358,6 @@ def zmat_process(self, zmat_output): [str(i - first_index + 1), List[1]], dtype=object ) self.bond_variables.append("R" + str(i - first_index)) - # self.angle_indices.append( - # [str(i - first_index + 1), List[1], List[2]] - # ) self.angle_indices = np.append(self.angle_indices, 0) self.angle_indices[-1] = np.array( [str(i - first_index + 1), List[1], List[2]], dtype=object @@ -424,6 +555,7 @@ def zmat_process(self, zmat_output): count += 1 new_walks = new_walks.reshape((-1, count)) + print(new_walks) new_walks = np.unique(new_walks, axis=0) del_list = np.array([]) diff --git a/examples/2_Methanol/molpro/CMA0A/manual_sym/main.py b/examples/2_Methanol/molpro/CMA0A/manual_sym/main.py index facfada..daf02e6 100644 --- a/examples/2_Methanol/molpro/CMA0A/manual_sym/main.py +++ b/examples/2_Methanol/molpro/CMA0A/manual_sym/main.py @@ -11,6 +11,8 @@ "energy_regex_b": r"\(T\) total energy\s+(\-\d+\.\d+)", "cart_insert_b": 9, "cart_insert_a": 9, + # "calc_b" : False, + # "gen_disps_b" : False, "man_proj": True, "coords": "Custom", "success_regex_b": r"Molpro calculation terminated", @@ -63,4 +65,5 @@ def normalize(mat): from concordantmodes.cma import ConcordantModes CMA_obj = ConcordantModes(options_obj, proj=Proj) +# CMA_obj.run() CMA_obj.run(sym_sort=sym_sort) diff --git a/examples/2_Methanol/psi4/CMA0A/gradients/main.py b/examples/2_Methanol/psi4/CMA0A/gradients/main.py index ebaf02c..d79318a 100644 --- a/examples/2_Methanol/psi4/CMA0A/gradients/main.py +++ b/examples/2_Methanol/psi4/CMA0A/gradients/main.py @@ -10,8 +10,8 @@ "gradient_regex_b": [r"Total Gradient", r"tstop"], "cart_insert_b": 7, "cart_insert_a": 7, - "gen_disps_b": False, - "calc_b": False, + # "gen_disps_b": False, + # "calc_b": False, "coords": "Delocalized", "success_regex_b": r"beer", "success_regex_a": r"beer",