From a4e6d18ab51c4cb5a0f8515ebd5d35b3ecd20e6f Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Sun, 31 Mar 2019 12:05:30 -0700 Subject: [PATCH 1/5] Handlers should not print newline by default --- geometric/internal.py | 93 +++++------- geometric/log.ini | 4 +- geometric/logJson.ini | 22 +++ geometric/molecule.py | 328 ++++++++++++++++++++++++++++++++++++++++-- geometric/nifty.py | 68 +++++---- geometric/optimize.py | 121 ++++++++-------- geometric/rotate.py | 28 ++-- geometric/run_json.py | 19 ++- 8 files changed, 508 insertions(+), 175 deletions(-) create mode 100644 geometric/logJson.ini diff --git a/geometric/internal.py b/geometric/internal.py index 5cfde95e..f8461a68 100644 --- a/geometric/internal.py +++ b/geometric/internal.py @@ -114,20 +114,6 @@ def d_nucross(a, b): return np.dot(d_unit_vector(a), d_ncross(ev, b)) ## End vector calculus functions -def logArray(mat, precision=3, fmt="f"): - fmt="%% .%i%s" % (precision, fmt) - if len(mat.shape) == 1: - for i in range(mat.shape[0]): - logger.info(fmt % mat[i]), - print - elif len(mat.shape) == 2: - for i in range(mat.shape[0]): - for j in range(mat.shape[1]): - logger.info(fmt % mat[i,j]), - print - else: - raise RuntimeError("One or two dimensional arrays only") - class CartesianX(object): def __init__(self, a, w=1.0): self.a = a @@ -143,7 +129,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = self.a == other.a if eq and self.w != other.w: - logger.warning("Warning: CartesianX same atoms, different weights (%.4f %.4f)" % (self.w, other.w)) + logger.warning("Warning: CartesianX same atoms, different weights (%.4f %.4f)\n" % (self.w, other.w)) return eq def __ne__(self, other): @@ -180,7 +166,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = self.a == other.a if eq and self.w != other.w: - logger.warning("Warning: CartesianY same atoms, different weights (%.4f %.4f)" % (self.w, other.w)) + logger.warning("Warning: CartesianY same atoms, different weights (%.4f %.4f)\n" % (self.w, other.w)) return eq def __ne__(self, other): @@ -217,7 +203,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = self.a == other.a if eq and self.w != other.w: - logger.warning("Warning: CartesianZ same atoms, different weights (%.4f %.4f)" % (self.w, other.w)) + logger.warning("Warning: CartesianZ same atoms, different weights (%.4f %.4f)\n" % (self.w, other.w)) return eq def __ne__(self, other): @@ -255,7 +241,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = set(self.a) == set(other.a) if eq and np.sum((self.w-other.w)**2) > 1e-6: - logger.warning("Warning: TranslationX same atoms, different weights") + logger.warning("Warning: TranslationX same atoms, different weights\n") eq = False return eq @@ -295,7 +281,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = set(self.a) == set(other.a) if eq and np.sum((self.w-other.w)**2) > 1e-6: - logger.warning("Warning: TranslationY same atoms, different weights") + logger.warning("Warning: TranslationY same atoms, different weights\n") eq = False return eq @@ -335,7 +321,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = set(self.a) == set(other.a) if eq and np.sum((self.w-other.w)**2) > 1e-6: - logger.warning("Warning: TranslationZ same atoms, different weights") + logger.warning("Warning: TranslationZ same atoms, different weights\n") eq = False return eq @@ -399,7 +385,7 @@ def __eq__(self, other): if type(self) is not type(other): return False eq = set(self.a) == set(other.a) if eq and np.sum((self.x0-other.x0)**2) > 1e-6: - logger.warning("Warning: Rotator same atoms, different reference positions") + logger.warning("Warning: Rotator same atoms, different reference positions\n") return eq def __repr__(self): @@ -1492,7 +1478,7 @@ def __eq__(self, other): if self.a == other.a: if {self.b, self.c, self.d} == {other.b, other.c, other.d}: if [self.b, self.c, self.d] != [other.b, other.c, other.d]: - logger.warning("Warning: OutOfPlane atoms are the same, ordering is different") + logger.warning("Warning: OutOfPlane atoms are the same, ordering is different\n") return True # if self.b == other.b: # if self.c == other.c: @@ -1622,7 +1608,7 @@ def wilsonB(self, xyz): WilsonB.append(Der[i].flatten()) self.stored_wilsonB[xhash] = np.array(WilsonB) if len(self.stored_wilsonB) > 1000 and not CacheWarning: - logger.warning("\x1b[91mWarning: more than 1000 B-matrices stored, memory leaks likely\x1b[0m") + logger.warning("\x1b[91mWarning: more than 1000 B-matrices stored, memory leaks likely\x1b[0m\n") CacheWarning = True ans = np.array(WilsonB) return ans @@ -1648,7 +1634,7 @@ def GInverse_SVD(self, xyz): U, S, VT = np.linalg.svd(G) time_svd = click() except np.linalg.LinAlgError: - logger.warning("\x1b[1;91m SVD fails, perturbing coordinates and trying again\x1b[0m") + logger.warning("\x1b[1;91m SVD fails, perturbing coordinates and trying again\x1b[0m\n") xyz = xyz + 1e-2*np.random.random(xyz.shape) loops += 1 if loops == 10: @@ -1693,7 +1679,7 @@ def checkFiniteDifferenceGrad(self, xyz): x2[i,j] -= h PMDiff = self.calcDiff(x1,x2) FiniteDifference[:,i,j] = PMDiff/(2*h) - logger.info("-=# Now checking first derivatives of internal coordinates w/r.t. Cartesians #=-") + logger.info("-=# Now checking first derivatives of internal coordinates w/r.t. Cartesians #=-\n") for i in range(Analytical.shape[0]): title = "%20s : %20s" % ("IC %i/%i" % (i+1, Analytical.shape[0]), self.Internals[i]) lines = [title] @@ -1710,9 +1696,9 @@ def checkFiniteDifferenceGrad(self, xyz): if maxerr < np.abs(error): maxerr = np.abs(error) if maxerr > 1e-5: - logger.info('\n'.join(lines)) - logger.info("%s : Max Error = %.5e" % (title, maxerr)) - logger.info("Finite-difference Finished") + logger.info('\n'.join(lines)+'\n') + logger.info("%s : Max Error = %.5e\n" % (title, maxerr)) + logger.info("Finite-difference Finished\n") return FiniteDifference def checkFiniteDifferenceHess(self, xyz): @@ -1721,7 +1707,7 @@ def checkFiniteDifferenceHess(self, xyz): FiniteDifference = np.zeros_like(Analytical) h = 1e-4 verbose = False - logger.info("-=# Now checking second derivatives of internal coordinates w/r.t. Cartesians #=-") + logger.info("-=# Now checking second derivatives of internal coordinates w/r.t. Cartesians #=-\n") for j in range(xyz.shape[0]): for m in range(3): for k in range(xyz.shape[0]): @@ -1746,7 +1732,7 @@ def checkFiniteDifferenceHess(self, xyz): for i in range(Analytical.shape[0]): title = "%20s : %20s" % ("IC %i/%i" % (i+1, Analytical.shape[0]), self.Internals[i]) lines = [title] - if verbose: logger.info(title) + if verbose: logger.info(title+'\n') maxerr = 0.0 numerr = 0 for j in range(Analytical.shape[1]): @@ -1761,14 +1747,14 @@ def checkFiniteDifferenceHess(self, xyz): if np.abs(error)>1e-5: numerr += 1 if (ana != 0.0 or fin != 0.0) and verbose: - logger.info(message) + logger.info(message+'\n') lines.append(message) if maxerr < np.abs(error): maxerr = np.abs(error) if maxerr > 1e-5 and not verbose: - logger.info('\n'.join(lines)) - logger.info("%s : Max Error = % 12.5e (%i above threshold)" % (title, maxerr, numerr)) - logger.info("Finite-difference Finished") + logger.info('\n'.join(lines)+'\n') + logger.info("%s : Max Error = % 12.5e (%i above threshold)\n" % (title, maxerr, numerr)) + logger.info("Finite-difference Finished\n") return FiniteDifference def calcGrad(self, xyz, gradx): @@ -1831,14 +1817,14 @@ def newCartesian(self, xyz, dQ, verbose=True): # Function to exit from loop def finish(microiter, rmsdt, ndqt, xyzsave, xyz_iter1): if ndqt > 1e-1: - if verbose: logger.info("Failed to obtain coordinates after %i microiterations (rmsd = %.3e |dQ| = %.3e)" % (microiter, rmsdt, ndqt)) + if verbose: logger.info("Failed to obtain coordinates after %i microiterations (rmsd = %.3e |dQ| = %.3e)\n" % (microiter, rmsdt, ndqt)) self.bork = True self.writeCache(xyz, dQ, xyz_iter1) return xyz_iter1.flatten() elif ndqt > 1e-3: - if verbose: logger.info("Approximate coordinates obtained after %i microiterations (rmsd = %.3e |dQ| = %.3e)" % (microiter, rmsdt, ndqt)) + if verbose: logger.info("Approximate coordinates obtained after %i microiterations (rmsd = %.3e |dQ| = %.3e)\n" % (microiter, rmsdt, ndqt)) else: - if verbose: logger.info("Cartesian coordinates obtained after %i microiterations (rmsd = %.3e |dQ| = %.3e)" % (microiter, rmsdt, ndqt)) + if verbose: logger.info("Cartesian coordinates obtained after %i microiterations (rmsd = %.3e |dQ| = %.3e)\n" % (microiter, rmsdt, ndqt)) self.writeCache(xyz, dQ, xyzsave) return xyzsave.flatten() fail_counter = 0 @@ -1858,19 +1844,19 @@ def finish(microiter, rmsdt, ndqt, xyzsave, xyz_iter1): ndq = np.linalg.norm(dQ1-dQ_actual) if len(ndqs) > 0: if ndq > ndqt: - if verbose: logger.info("Iter: %i Err-dQ (Best) = %.5e (%.5e) RMSD: %.5e Damp: %.5e (Bad)" % (microiter, ndq, ndqt, rmsd, damp)) + if verbose: logger.info("Iter: %i Err-dQ (Best) = %.5e (%.5e) RMSD: %.5e Damp: %.5e (Bad)\n" % (microiter, ndq, ndqt, rmsd, damp)) damp /= 2 fail_counter += 1 # xyz2 = xyz1.copy() else: - if verbose: logger.info("Iter: %i Err-dQ (Best) = %.5e (%.5e) RMSD: %.5e Damp: %.5e (Good)" % (microiter, ndq, ndqt, rmsd, damp)) + if verbose: logger.info("Iter: %i Err-dQ (Best) = %.5e (%.5e) RMSD: %.5e Damp: %.5e (Good)\n" % (microiter, ndq, ndqt, rmsd, damp)) fail_counter = 0 damp = min(damp*1.2, 1.0) rmsdt = rmsd ndqt = ndq xyzsave = xyz2.copy() else: - if verbose: logger.info("Iter: %i Err-dQ = %.5e RMSD: %.5e Damp: %.5e" % (microiter, ndq, rmsd, damp)) + if verbose: logger.info("Iter: %i Err-dQ = %.5e RMSD: %.5e Damp: %.5e\n" % (microiter, ndq, rmsd, damp)) rmsdt = rmsd ndqt = ndq ndqs.append(ndq) @@ -2258,14 +2244,14 @@ def update(self, other): else: i.inactive = 0 if i.inactive == 1: - logger.info("Deleting:", i) + logger.info("Deleting:" + str(i) + "\n") self.Internals.remove(i) Changed = True else: i.inactive = 0 for i in other.Internals: if i not in self.Internals: - logger.info("Adding: ", i) + logger.info("Adding: " + str(i) + "\n") self.Internals.append(i) Changed = True return Changed @@ -2274,7 +2260,7 @@ def join(self, other): Changed = False for i in other.Internals: if i not in self.Internals: - logger.info("Adding: ", i) + logger.info("Adding: " + str(i) + "\n") self.Internals.append(i) Changed = True return Changed @@ -2356,13 +2342,13 @@ def getRotatorDots(self): def printRotations(self, xyz): rotNorms = self.getRotatorNorms() if len(rotNorms) > 0: - logger.info("Rotator Norms: " + " ".join(["% .4f" % i for i in rotNorms])) + logger.info("Rotator Norms: " + " ".join(["% .4f" % i for i in rotNorms]) + "\n") rotDots = self.getRotatorDots() if len(rotDots) > 0 and np.max(rotDots) > 1e-5: - logger.info("Rotator Dots : " + " ".join(["% .4f" % i for i in rotDots])) + logger.info("Rotator Dots : " + " ".join(["% .4f" % i for i in rotDots]) + "\n") linAngs = [ic.value(xyz) for ic in self.Internals if type(ic) is LinearAngle] if len(linAngs) > 0: - logger.info("Linear Angles: " + " ".join(["% .4f" % i for i in linAngs])) + logger.info("Linear Angles: " + " ".join(["% .4f" % i for i in linAngs]) + "\n") def derivatives(self, xyz): self.calculate(xyz) @@ -2426,7 +2412,7 @@ def addConstraint(self, cPrim, cVal=None, xyz=None): if cPrim in self.cPrims: iPrim = self.cPrims.index(cPrim) if np.abs(cVal - self.cVals[iPrim]) > 1e-6: - logger.info("Updating constraint value to %.4e" % cVal) + logger.info("Updating constraint value to %.4e\n" % cVal) self.cVals[iPrim] = cVal else: if cPrim not in self.Internals: @@ -2498,11 +2484,10 @@ def printConstraints(self, xyz, thre=1e-5): if np.abs(diff*factor) > thre: out_lines.append("%-30s % 10.5f % 10.5f % 10.5f" % (str(c), current*factor, reference*factor, diff*factor)) if len(out_lines) > 0: - logger.info(header) - logger.info('\n'.join(out_lines)) + logger.info(header + "\n") + logger.info('\n'.join(out_lines) + "\n") # if type(c) in [RotationA, RotationB, RotationC]: # print c, c.value(xyz) - # logArray(c.x0) def getConstraintTargetVals(self): nc = len(self.cPrims) @@ -2764,7 +2749,7 @@ def applyConstraints(self, xyz): if niter > 1 and np.linalg.norm(dQ) > np.linalg.norm(dQ0): xyz1 = xyzs[np.argmin(ndqs)] if not self.enforce_fail_printed: - logger.warning("Warning: Failed to enforce exact constraint satisfaction. Please remove possible redundant constraints. See below:") + logger.warning("Warning: Failed to enforce exact constraint satisfaction. Please remove possible redundant constraints. See below:\n") self.printConstraints(xyz1, thre=0.0) self.enforce_fail_printed = True return xyz1 @@ -2789,7 +2774,7 @@ def newCartesian_withConstraint(self, xyz, dQ, thre=0.1, verbose=False): if constraintSmall: xyz2 = self.applyConstraints(xyz2) if not self.enforced: - logger.info("<<< Enforcing constraint satisfaction >>>") + logger.info("<<< Enforcing constraint satisfaction >>>\n") self.enforced = True else: self.enforced = False @@ -3032,7 +3017,7 @@ def ov(vi, vj): ui = U[:, ic] Unorms[ic] = np.sqrt(ov(ui, ui)) if Unorms[ic]/Vnorms[ic] < 0.1: - logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original" % (ic, Unorms[ic]/Vnorms[ic])) + logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original\n" % (ic, Unorms[ic]/Vnorms[ic])) V0 = V.copy() # Project out newest U column from all remaining V columns. for jc in range(ic+1, nv): @@ -3144,7 +3129,7 @@ def ov(vi, vj): ui = U[:, ic] Unorms[ic] = np.sqrt(ov(ui, ui)) if Unorms[ic]/Vnorms[ic] < 0.1: - logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original" % (ic-6, Unorms[ic]/Vnorms[ic])) + logger.warning("Constraint %i is almost redundant; after projection norm is %.3f of original\n" % (ic-6, Unorms[ic]/Vnorms[ic])) V0 = V.copy() # Project out newest U column from all remaining V columns. for jc in range(ic+1, nv): diff --git a/geometric/log.ini b/geometric/log.ini index b27c9fe5..0fca0beb 100644 --- a/geometric/log.ini +++ b/geometric/log.ini @@ -12,13 +12,13 @@ level=INFO handlers=stream_handler, file_handler [handler_stream_handler] -class=StreamHandler +class=geometric.nifty.RawStreamHandler level=INFO formatter=formatter args=(sys.stderr,) [handler_file_handler] -class=FileHandler +class=geometric.nifty.RawFileHandler level=INFO formatter=formatter args=('%(logfilename)s',) diff --git a/geometric/logJson.ini b/geometric/logJson.ini new file mode 100644 index 00000000..626ee17f --- /dev/null +++ b/geometric/logJson.ini @@ -0,0 +1,22 @@ +[loggers] +keys=root + +[handlers] +keys=stream_handler + +[formatters] +keys=formatter + +[logger_root] +level=INFO +handlers=stream_handler + +[handler_stream_handler] +class=geometric.nifty.RawStreamHandler +level=INFO +formatter=formatter +args=(sys.stderr,) + +[formatter_formatter] +format=%(message)s +#format=%(asctime)s %(name)-12s %(levelname)-8s %(message)s diff --git a/geometric/molecule.py b/geometric/molecule.py index b2f67259..2bbc0a9b 100644 --- a/geometric/molecule.py +++ b/geometric/molecule.py @@ -37,7 +37,7 @@ # | Chemical file format conversion module |# # | |# # | Lee-Ping Wang (leeping@ucdavis.edu) |# -# | Last updated June 6, 2018 |# +# | Last updated March 31, 2019 |# # | |# # | This code is part of geomeTRIC and is covered under the |# # | geomeTRIC copyright notice and MIT license. |# @@ -188,9 +188,11 @@ if "forcebalance" in __name__: # If this module is part of ForceBalance, use the package level logger from .output import * + package = "ForceBalance" elif "geometric" in __name__: # This ensures logging behavior is consistent with the rest of geomeTRIC from .nifty import logger + package = "geomeTRIC" else: # Previous default behavior if FB package level loggers could not be imported from logging import * @@ -213,6 +215,7 @@ def emit(self, record): logger.setLevel(INFO) handler = RawStreamHandler() logger.addHandler(handler) + package = "molecule.py" module_name = __name__.replace('.molecule','') @@ -1013,6 +1016,76 @@ def AtomContact(xyz, pairs, box=None, displace=False): else: return dr +#===================================# +#| Rotation subroutine 2018-10-28 |# +#| Copied from geomeTRIC.rotate |# +#===================================# + +def form_rot(q): + """ + Given a quaternion p, form a rotation matrix from it. + + Parameters + ---------- + q : numpy.ndarray + 1D array with 3 elements representing the rotation quaterion. + Elements of quaternion are : [cos(a/2), sin(a/2)*axis[0..2]] + + Returns + ------- + numpy.array + 3x3 rotation matrix + """ + assert q.ndim == 1 + assert q.shape[0] == 4 + # Take the "complex conjugate" + qc = np.zeros_like(q) + qc[0] = q[0] + qc[1] = -q[1] + qc[2] = -q[2] + qc[3] = -q[3] + # Form al_q and al_qc matrices + al_q = np.array([[ q[0], -q[1], -q[2], -q[3]], + [ q[1], q[0], -q[3], q[2]], + [ q[2], q[3], q[0], -q[1]], + [ q[3], -q[2], q[1], q[0]]]) + ar_qc = np.array([[ qc[0], -qc[1], -qc[2], -qc[3]], + [ qc[1], qc[0], qc[3], -qc[2]], + [ qc[2], -qc[3], qc[0], qc[1]], + [ qc[3], qc[2], -qc[1], qc[0]]]) + # Multiply matrices + R4 = np.dot(al_q,ar_qc) + return R4[1:, 1:] + +def axis_angle(axis, angle): + """ + Given a rotation axis and angle, return the corresponding + 3x3 rotation matrix, which will rotate a (Nx3) array of + xyz coordinates as x0_rot = np.dot(R, x0.T).T + + Parameters + ---------- + axis : numpy.ndarray + 1D array with 3 elements representing the rotation axis + angle : float + The angle of the rotation + + Returns + ------- + numpy.array + 3x3 rotation matrix + """ + assert axis.ndim == 1 + assert axis.shape[0] == 3 + axis /= np.linalg.norm(axis) + # Make quaternion + ct2 = np.cos(angle/2) + st2 = np.sin(angle/2) + q = np.array([ct2, st2*axis[0], st2*axis[1], st2*axis[2]]) + # Form rotation matrix + R = form_rot(q) + return R + class Molecule(object): """ Lee-Ping's general file format conversion class. @@ -1540,10 +1613,10 @@ def write(self,fnm=None,ftype=None,append=False,selection=None,**kwargs): ftype = os.path.splitext(fnm)[1][1:] ## Fill in comments. if 'comms' not in self.Data: - self.comms = ['Generated by ForceBalance from %s: Frame %i of %i' % (fnm, i+1, self.ns) for i in range(self.ns)] + self.comms = ['Generated by %s from %s: Frame %i of %i' % (package, fnm, i+1, self.ns) for i in range(self.ns)] if 'xyzs' in self.Data and len(self.comms) < len(self.xyzs): for i in range(len(self.comms), len(self.xyzs)): - self.comms.append("Frame %i: generated by ForceBalance" % i) + self.comms.append("Frame %i: generated by %s" % (i, package)) ## I needed to add in this line because the DCD writer requires the file name, ## but the other methods don't. self.fout = fnm @@ -1645,7 +1718,7 @@ def load_frames(self, fnm, ftype=None, **kwargs): self.Data[key] = val ## Create a list of comment lines if we don't already have them from reading the file. if 'comms' not in self.Data: - self.comms = ['Generated by ForceBalance from %s: Frame %i of %i' % (fnm, i+1, self.ns) for i in range(self.ns)] + self.comms = ['Generated by %s from %s: Frame %i of %i' % (package, fnm, i+1, self.ns) for i in range(self.ns)] else: self.comms = [i.expandtabs() for i in self.comms] @@ -2105,6 +2178,237 @@ def distance_displacement(self): dxij.append(dxij_i) return AtomIterator, drij, dxij + def rotate_bond(self, frame, aj, ak, increment=15): + """ + Return a new Molecule object containing the selected frame + plus a number of frames where the selected dihedral angle is rotated + in steps of 'increment' given in degrees. + + This function is designed to be called by Molecule.rotate_check_clash(). + + Parameters + ---------- + frame : int + Structure number of the current Molecule to be rotated + aj, ak : int + Atom numbers of the bond to be rotated + increment : float + Degrees of the rotation increment + + Returns + ------- + Molecule + New Molecule object containing the rotated structures + """ + # Select the single frame containing the structure to be rotated. + M = self[frame] + + # Delete the connection between atoms to be rotated in order + # to determine the rotation fragments. + delBonds = [] + for ibond, bond in enumerate(M.bonds): + if bond == (aj, ak) or bond == (ak, aj): + delBonds.append(bond) + if len(delBonds) > 1: + raise RuntimeError('Expected only one bond to be deleted') + if len(M.molecules) != 1: + raise RuntimeError('Expected a single molecule') + M.bonds.remove(delBonds[0]) + M.top_settings['read_bonds']=True + M.build_topology(force_bonds=False) + if len(M.molecules) != 2: + raise RuntimeError('Expected two molecules after removing a bond') + + # gAtoms contains the set of atoms to be rotated + # oAtoms contains the "other" atoms + if len(M.molecules[0].L()) < len(M.molecules[1].L()): + gAtoms = M.molecules[0].L() + oAtoms = M.molecules[1].L() + else: + gAtoms = M.molecules[1].L() + oAtoms = M.molecules[0].L() + + # atom2 is the "reference atom" on the group being rotated + # and will be moved to the origin during the rotation operation + if aj in gAtoms: + atom1 = ak + atom2 = aj + else: + atom1 = aj + atom2 = ak + M.bonds.append(delBonds[0]) + + # Rotation axis + axis = M.xyzs[0][atom2] - M.xyzs[0][atom1] + + # Move the "reference atom" to the origin + x0 = M.xyzs[0][gAtoms] + x0_ref = M.xyzs[0][atom2] + x0 -= x0_ref + + # Create grid in rotation angle + # and the rotated structures + for thetaDeg in np.arange(increment, 360, increment): + theta = np.pi * thetaDeg / 180 + # Make quaternion + R = axis_angle(axis, theta) + # Get rotated coordinates + x0_rot = np.dot(R, x0.T).T + # Copy old coordinates to new + xnew = M.xyzs[0].copy() + # Write rotated positions into new coordinates; + # shift back to original positions + xnew[gAtoms] = x0_rot + x0_ref + # Append new coordinates to our Molecule object + M.xyzs.append(xnew) + M.comms.append("Rotated by %.2f degrees" % increment) + return M, (gAtoms, oAtoms) + + def find_clashes(self, thre=0.0, pbc=True, groups=None): + """ + Obtain a list of atoms that 'clash' (i.e. are more than + 3 bonds apart and are closer than the provided threshold.) + + Parameters + ---------- + thre : float + Create a sorted-list of all non-bonded atom pairs + with distance below this threshold + pbc : bool + Whether to use PBC when computing interatomic distances + filt : tuple (group1, group2) + If not None, only compute distances between atoms in 'group1' + and atoms in 'group2'. For example this could be [gAtoms, oAtoms] + from rotate_bond + + Returns + ------- + minPair_frames : list of tuples + The closest pair of non-bonded atoms in each frame + minDist_frames : list of tuples + The distance between the closest pair of non-bonded atoms in each frame + clashPairs_frames : list of lists + Pairs of non-bonded atoms with distance below the threshold in each frame + clashDists_frames : list of lists + Distances between pairs of non-bonded atoms below the threshold in each frame + """ + if groups is not None: + AtomIterator = np.ascontiguousarray([[min(g), max(g)] for g in itertools.product(groups[0], groups[1])]) + else: + AtomIterator = np.ascontiguousarray(np.vstack((np.fromiter(itertools.chain(*[[i]*(self.na-i-1) for i in range(self.na)]),dtype=np.int32), + np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) + ang13 = [(min(a[0], a[2]), max(a[0], a[2])) for a in self.find_angles()] + dih14 = [(min(d[0], d[3]), max(d[0], d[3])) for d in self.find_dihedrals()] + bondedPairs = np.where([tuple(aPair) in (self.bonds+ang13+dih14) for aPair in AtomIterator])[0] + AtomIterator_nb = np.delete(AtomIterator, bondedPairs, axis=0) + + minPair_frames = [] + minDist_frames = [] + clashPairs_frames = [] + clashDists_frames = [] + for frame in range(len(self)): + if hasattr(self, 'boxes') and pbc: + box=np.array([self.boxes[frame].a, self.boxes[frame].b, self.boxes[frame].c]) + drij = AtomContact(self.xyzs[frame],AtomIterator_nb,box=box) + else: + drij = AtomContact(self.xyzs[frame],AtomIterator_nb) + clashPairIdx = np.where(drij < thre)[0] + clashPairs = AtomIterator_nb[clashPairIdx] + clashDists = drij[clashPairIdx] + sorter = np.argsort(clashDists) + clashPairs = clashPairs[sorter] + clashDists = clashDists[sorter] + minIdx = np.argmin(drij) + minPair = AtomIterator_nb[minIdx] + minDist = drij[minIdx] + minPair_frames.append(minPair) + minDist_frames.append(minDist) + clashPairs_frames.append(clashPairs.copy()) + clashDists_frames.append(clashDists.copy()) + return minPair_frames, minDist_frames, clashPairs_frames, clashDists_frames + + def rotate_check_clash(self, frame, rotate_index, thresh_hyd=1.4, thresh_hvy=1.8, printLevel=1): + """ + Return a new Molecule object containing the selected frame + plus a number of frames where the selected dihedral angle is rotated + in steps of 'increment' given in degrees. Additionally, check for + if pairs of non-bonded atoms "clash" i.e. approach below the specified + thresholds. + + Parameters + ---------- + frame : int + Structure number of the current Molecule to be rotated + (if only a single frame, this number should be 0) + rotate_index : tuple + 4 atom indices (ai, aj, ak, al) representing the dihedral angle to be rotated. + The first and last indices are used to measure the dihedral angle. + thresh_hyd : float + Clash threshold for hydrogen atoms. Reasonable values are in between 1.3 and 2.0. + thresh_hvy : float + Clash threshold for heavy atoms. Reasonable values are in between 1.7 and 2.5. + printLevel: int + Sets the amount of printout (larger = more printout) + + Returns + ------- + Molecule + New Molecule object containing the rotated structures + Success : bool + True if no clashes + """ + if not hasattr(self, 'atomname'): + raise RuntimeError("Please add atom names before calling rotate_check_clash().") + ai, aj, ak, al = rotate_index + Success = False + # Create grid of rotated structures + M_rot_H, frags = self.rotate_bond(frame, aj, ak, 15) + phis = M_rot_H.measure_dihedrals(*rotate_index) + for i in range(len(M_rot_H)): + M_rot_H.comms[i] = ('Rigid scan: atomname %s, serial %s, dihedral %.3f' + % ('-'.join([self.atomname[i] for i in rotate_index]), + '-'.join(["%i" % (i+1) for i in rotate_index]), phis[i])) + heavyIdx = [i for i in range(self.na) if self.elem[i] != 'H'] + heavy_frags = [[],[]] + for iHeavy, iAll in enumerate(heavyIdx): + if iAll in frags[0]: + heavy_frags[0].append(iHeavy) + if iAll in frags[1]: + heavy_frags[1].append(iHeavy) + M_rot_C = M_rot_H.atom_select(heavyIdx) + minPair_H_frames, minDist_H_frames, clashPairs_H_frames, clashDists_H_frames = M_rot_H.find_clashes(thre=thresh_hyd, groups=frags) + minPair_C_frames, minDist_C_frames, clashPairs_C_frames, clashDists_C_frames = M_rot_C.find_clashes(thre=thresh_hvy, groups=heavy_frags) + # Get the following information: (1) Whether a clash exists, (2) the frame with the smallest distance, + # (3) the pair of atoms with the smallest distance, (4) the smallest distance + haveClash_H = any([len(c) > 0 for c in clashPairs_H_frames]) + minFrame_H = np.argmin(minDist_H_frames) + minAtoms_H = minPair_H_frames[minFrame_H] + minDist_H = minDist_H_frames[minFrame_H] + haveClash_C = any([len(c) > 0 for c in clashPairs_C_frames]) + minFrame_C = np.argmin(minDist_C_frames) + minAtoms_C = minPair_C_frames[minFrame_C] + minDist_C = minDist_C_frames[minFrame_C] + if not (haveClash_H or haveClash_C): + if printLevel >= 1: + print("\n \x1b[1;92mSuccess - no clashes. Thresh(H, Hvy) = (%.2f, %.2f)\x1b[0m" % (thresh_hyd, thresh_hvy)) + mini = M_rot_H.atomname[minAtoms_H[0]] + minj = M_rot_H.atomname[minAtoms_H[1]] + print(" Closest (Hyd) : rot-frame %i atoms %s-%s %.2f" % (minFrame_H, mini, minj, minDist_H)) + mini = M_rot_C.atomname[minAtoms_C[0]] + minj = M_rot_C.atomname[minAtoms_C[1]] + print(" Closest (Hvy) : rot-frame %i atoms %s-%s %.2f" % (minFrame_C, mini, minj, minDist_C)) + Success = True + else: + if haveClash_H: + mini = M_rot_H.atomname[minAtoms_H[0]] + minj = M_rot_H.atomname[minAtoms_H[1]] + if printLevel >= 2: print(" Clash (Hyd) : rot-frame %i atoms %s-%s %.2f" % (minFrame_H, mini, minj, minDist_H)) + if haveClash_C: + mini = M_rot_C.atomname[minAtoms_C[0]] + minj = M_rot_C.atomname[minAtoms_C[1]] + if printLevel >= 2: print(" Clash (Hvy) : rot-frame %i atoms %s-%s %.2f" % (minFrame_C, mini, minj, minDist_C)) + return M_rot_H, Success + def find_angles(self): """ Return a list of 3-tuples corresponding to all of the @@ -3291,7 +3595,8 @@ def read_pdb(self, fnm, **kwargs): conect_B_list.append(int(line_rest[:5]) - 1) line_rest = line_rest[5:] for conect_B in conect_B_list: - bonds.append([conect_A, conect_B]) + bond = (min((conect_A, conect_B)), max((conect_A, conect_B))) + bonds.append(bond) Answer={"xyzs":XYZList, "chain":list(ChainID), "altloc":list(AltLoc), "icode":list(ICode), "atomname":[str(i) for i in AtomNames], "resid":list(ResidueID), "resname":list(ResidueNames), @@ -3881,7 +4186,7 @@ def write_molproq(self, selection, **kwargs): def write_mdcrd(self, selection, **kwargs): self.require('xyzs') # In mdcrd files, there is only one comment line - out = ['mdcrd file generated using ForceBalance'] + out = ['mdcrd file generated using %s' % package] for I in selection: xyz = self.xyzs[I] out += [''.join(["%8.3f" % i for i in g]) for g in grouper(10, list(xyz.flatten()))] @@ -3955,7 +4260,6 @@ def write_gro(self, selection, **kwargs): xyzwrite = xyz.copy() xyzwrite /= 10.0 # GROMACS uses nanometers out.append(self.comms[I]) - #out.append("Generated by ForceBalance from %s" % self.fnm) out.append("%5i" % self.na) for an, line in enumerate(xyzwrite): out.append(format_gro_coord(self.resid[an],self.resname[an],atomname[an],an+1,xyzwrite[an])) @@ -4056,7 +4360,7 @@ def write_pdb(self, selection, **kwargs): out = [] # Create the PDB header. - out.append("REMARK 1 CREATED WITH FORCEBALANCE %s" % (str(date.today()))) + out.append("REMARK 1 CREATED WITH %s %s" % (package.upper(), str(date.today()))) if 'boxes' in self.Data: a = self.boxes[0].a b = self.boxes[0].b @@ -4080,10 +4384,14 @@ def write_pdb(self, selection, **kwargs): resId = resIds[i] coords = self.xyzs[sn][i] symbol = self.elem[i] + if hasattr(self, 'partial_charge'): + bfactor = self.partial_charge[i] + else: + bfactor = 0.0 atomIndices[i] = atomIndex - line = "%s%5d %-4s %3s %s%4s %s%s%s 1.00 0.00 %2s " % ( + line = "%s%5d %-4s %3s %s%4s %s%s%s %5.2f 0.00 %2s " % ( recordName, atomIndex%100000, atomName, resName, chainName, resId, _format_83(coords[0]), - _format_83(coords[1]), _format_83(coords[2]), symbol) + _format_83(coords[1]), _format_83(coords[2]), bfactor, symbol) assert len(line) == 80, 'Fixed width overflow detected' out.append(line) atomIndex += 1 diff --git a/geometric/nifty.py b/geometric/nifty.py index d777a396..f3b2e4e1 100644 --- a/geometric/nifty.py +++ b/geometric/nifty.py @@ -50,19 +50,16 @@ if "forcebalance" in __name__: # If this module is part of ForceBalance, use the package level logger from .output import * -elif "geometric" in __name__: - # This ensures logging behavior is consistent with the rest of geomeTRIC - from logging import * - logger = getLogger(__name__) - logger.setLevel(INFO) + package="ForceBalance" else: - # Previous default behavior if FB package level loggers could not be imported from logging import * + # Define two handlers that don't print newline characters at the end of each line class RawStreamHandler(StreamHandler): - """Exactly like output.StreamHandler except it does no extra formatting - before sending logging messages to the stream. This is more compatible with - how output has been displayed in ForceBalance. Default stream has also been - changed from stderr to stdout""" + """ + Exactly like StreamHandler, except no newline character is printed at the end of each message. + This is done in order to ensure functions in molecule.py and nifty.py work consistently + across multiple packages. + """ def __init__(self, stream = sys.stdout): super(RawStreamHandler, self).__init__(stream) @@ -70,14 +67,31 @@ def emit(self, record): message = record.getMessage() self.stream.write(message) self.flush() - # logger=getLogger() - # logger.handlers = [RawStreamHandler(sys.stdout)] - # LPW: Daniel Smith suggested these changes to improve logger behavior - logger = getLogger("NiftyLogger") - logger.setLevel(INFO) - handler = RawStreamHandler() - logger.addHandler(handler) + class RawFileHandler(FileHandler): + """ + Exactly like FileHandler, except no newline character is printed at the end of each message. + This is done in order to ensure functions in molecule.py and nifty.py work consistently + across multiple packages. + """ + def __init__(self, *args, **kwargs): + super(RawFileHandler, self).__init__(*args, **kwargs) + def emit(self, record): + if self.stream is None: + self.stream = self._open() + RawStreamHandler.emit(self, record) + + if "geometric" in __name__: + # This ensures logging behavior is consistent with the rest of geomeTRIC + logger = getLogger(__name__) + logger.setLevel(INFO) + package="geomeTRIC" + else: + logger = getLogger("NiftyLogger") + logger.setLevel(INFO) + handler = RawStreamHandler() + logger.addHandler(handler) + package="nifty.py" try: import bz2 @@ -93,7 +107,6 @@ def emit(self, record): logger.warning("gzip module import failed (used in compressing or decompressing pickle files)\n") HaveGZ = False - ## Boltzmann constant kb = 0.0083144100163 ## Q-Chem to GMX unit conversion for energy @@ -152,6 +165,7 @@ def pvec1d(vec1d, precision=1, format="e", loglevel=INFO): v2a = np.array(vec1d) for i in range(v2a.shape[0]): logger.log(loglevel, "%% .%i%s " % (precision, format) % v2a[i]) + logger.log(loglevel, '\n') def astr(vec1d, precision=4): """ Write an array to a string so we can use it to key a dictionary. """ @@ -164,7 +178,9 @@ def pmat2d(mat2d, precision=1, format="e", loglevel=INFO): """ m2a = np.array(mat2d) for i in range(m2a.shape[0]): - logger.log(loglevel, ' '.join(["%% .%i%s " % (precision, format) % m2a[i][j] for j in range(m2a[i].shape[0])])) + for j in range(m2a.shape[1]): + logger.log(loglevel, "%% .%i%s " % (precision, format) % m2a[i][j]) + logger.log(loglevel, '\n') def grouper(iterable, n): """Collect data into fixed-length chunks or blocks""" @@ -504,16 +520,16 @@ def monotonic_decreasing(arr, start=None, end=None, verbose=False): end = len(arr) - 1 a0 = arr[start] idx = [start] - if verbose: logger.info("Starting @ %i : %.6f" % (start, arr[start])) + if verbose: logger.info("Starting @ %i : %.6f\n" % (start, arr[start])) if end > start: i = start+1 while i < end: if arr[i] < a0: a0 = arr[i] idx.append(i) - if verbose: logger.info("Including %i : %.6f" % (i, arr[i])) + if verbose: logger.info("Including %i : %.6f\n" % (i, arr[i])) else: - if verbose: logger.info("Excluding %i : %.6f" % (i, arr[i])) + if verbose: logger.info("Excluding %i : %.6f\n" % (i, arr[i])) i += 1 if end < start: i = start-1 @@ -521,9 +537,9 @@ def monotonic_decreasing(arr, start=None, end=None, verbose=False): if arr[i] < a0: a0 = arr[i] idx.append(i) - if verbose: logger.info("Including %i : %.6f" % (i, arr[i])) + if verbose: logger.info("Including %i : %.6f\n" % (i, arr[i])) else: - if verbose: logger.info("Excluding %i : %.6f" % (i, arr[i])) + if verbose: logger.info("Excluding %i : %.6f\n" % (i, arr[i])) i -= 1 return np.array(idx) @@ -835,7 +851,7 @@ def getWQIds(): global WQIDS return WQIDS -def createWorkQueue(wq_port, debug=True, name='geomeTRIC'): +def createWorkQueue(wq_port, debug=True, name=package): global WORK_QUEUE if debug: work_queue.set_debug_flag('all') @@ -1083,7 +1099,7 @@ def listfiles(fnms=None, ext=None, err=False, dnm=None): raise RuntimeError answer = [fnms] elif fnms is not None: - logger.info(fnms) + logger.info(str(fnms)) logger.error('First argument to listfiles must be a list, a string, or None') raise RuntimeError if answer == [] and ext is not None: diff --git a/geometric/optimize.py b/geometric/optimize.py index d7ea895b..69c83e35 100644 --- a/geometric/optimize.py +++ b/geometric/optimize.py @@ -22,7 +22,6 @@ from .errors import EngineError, GeomOptNotConvergedError from enum import Enum - def RebuildHessian(IC, H0, coord_seq, grad_seq, params): """ Rebuild the Hessian after making a change to the internal coordinate system. @@ -57,7 +56,7 @@ def RebuildHessian(IC, H0, coord_seq, grad_seq, params): history += 1 if history < 1: return H0.copy() - logger.info("Rebuilding Hessian using %i gradients" % history) + logger.info("Rebuilding Hessian using %i gradients\n" % history) y_seq = [IC.calculate(i) for i in coord_seq[-history-1:]] g_seq = [IC.calcGrad(i, j) for i, j in zip(coord_seq[-history-1:],grad_seq[-history-1:])] Yprev = y_seq[0] @@ -77,7 +76,7 @@ def RebuildHessian(IC, H0, coord_seq, grad_seq, params): Hstor = H.copy() H += Mat1-Mat2 if np.min(np.linalg.eigh(H)[0]) < params.epsilon and params.reset: - logger.info("Eigenvalues below %.4e (%.4e) - returning guess" % (params.epsilon, np.min(np.linalg.eigh(H)[0]))) + logger.info("Eigenvalues below %.4e (%.4e) - returning guess\n" % (params.epsilon, np.min(np.linalg.eigh(H)[0]))) return H0.copy() return H @@ -226,7 +225,7 @@ def brent_wiki(f, a, b, rel, cvg=0.1, obj=None, verbose=False): # Convergence failure - interval becomes # smaller than threshold if np.abs(b-a) < epsilon: - if verbose: logger.info("returning because interval is too small") + if verbose: logger.info("returning because interval is too small\n") if obj is not None: obj.brentFailed = True return s # Exit before converging when @@ -245,11 +244,6 @@ def brent_wiki(f, a, b, rel, cvg=0.1, obj=None, verbose=False): a, b = b, a fa, fb = fb, fa -def ftest(x): - answer = (x+3)*(x-1)**2 - logger.info("(x, y) = ", x, answer) - return answer - def OneDScan(init, final, steps): """ Return a list of N equally spaced values between initial and final. @@ -327,7 +321,7 @@ def ParseConstraints(molecule, constraints_string): # This is a list-of-lists. The intention is to create a multidimensional grid # of constraint values if necessary. if len(line) == 0: continue - logger.info(line) + logger.info(line+'\n') if line.startswith("$"): mode = line.replace("$","") else: @@ -455,7 +449,7 @@ def ParseConstraints(molecule, constraints_string): # Get the angle theta1 = float(s[5]) * np.pi / 180 if np.abs(theta1) > np.pi * 0.9: - logger.info("Large rotation: Your constraint may not work") + logger.info("Large rotation: Your constraint may not work\n") if mode == "set": c = np.cos(theta1/2.0) s = np.sin(theta1/2.0) @@ -468,7 +462,7 @@ def ParseConstraints(molecule, constraints_string): elif mode == "scan": theta2 = float(s[6]) * np.pi / 180 if np.abs(theta2) > np.pi * 0.9: - logger.info("Large rotation: Your constraint may not work") + logger.info("Large rotation: Your constraint may not work\n") steps = int(s[7]) # To alleviate future confusion: # There is one group of three constraints that we are going to scan over in one dimension. @@ -531,11 +525,11 @@ def get_delta_prime_trm(v, X, G, H, IC, verbose=False): HT[i, i] = 0.0 if verbose >= 2: seig = sorted(np.linalg.eig(HT)[0]) - logger.info("sorted(eig) : % .5e % .5e % .5e ... % .5e % .5e % .5e" % (seig[0], seig[1], seig[2], seig[-3], seig[-2], seig[-1])) + logger.info("sorted(eig) : % .5e % .5e % .5e ... % .5e % .5e % .5e\n" % (seig[0], seig[1], seig[2], seig[-3], seig[-2], seig[-1])) try: Hi = invert_svd(HT) except: - logger.info("\x1b[1;91mSVD Error - increasing v by 0.001 and trying again\x1b[0m") + logger.info("\x1b[1;91mSVD Error - increasing v by 0.001 and trying again\x1b[0m\n") return get_delta_prime_trm(v+0.001, X, G, H, IC) dyc = flat(-1 * np.dot(Hi,col(GC))) dy = dyc[:len(G)] @@ -708,7 +702,7 @@ def trust_step(target, v0, X, G, H, IC, rfo, verbose=False): dy, sol, dy_prime, = get_delta_prime(v, X, G, H, IC, rfo, verbose) ndy = np.linalg.norm(dy) # This is "too" verbose; will enable integer print level later - if verbose >= 2: logger.info("v = %.5f dy -> target = %.5f -> %.5f" % (v, ndy, target)) + if verbose >= 2: logger.info("v = %.5f dy -> target = %.5f -> %.5f\n" % (v, ndy, target)) if np.abs((ndy-target)/target) < 0.001: return dy, sol # With Lagrange multipliers it may be impossible to go under a target step size @@ -722,10 +716,10 @@ def trust_step(target, v0, X, G, H, IC, rfo, verbose=False): m_sol = sol # Break out of infinite oscillation loops if niter%100 == 99: - logger.info("trust_step hit niter = 100, randomizing") + logger.info("trust_step hit niter = 100, randomizing\n") v += np.random.random() * niter / 100 if niter%1000 == 999: - logger.info("trust_step hit niter = 1000, giving up") + logger.info("trust_step hit niter = 1000, giving up\n") return m_dy, m_sol class Froot(object): @@ -789,7 +783,7 @@ def evaluate(self, trial): if self.stored_val is None or cnorm > self.stored_val: self.stored_arg = trial self.stored_val = cnorm - if self.params.verbose: logger.info("dy(i): %.4f dy(c) -> target: %.4f -> %.4f%s" % (trial, cnorm, self.target, " (done)" if self.from_above else "")) + if self.params.verbose: logger.info("dy(i): %.4f dy(c) -> target: %.4f -> %.4f%s\n" % (trial, cnorm, self.target, " (done)" if self.from_above else "")) return cnorm-self.target class OptParams(object): @@ -948,9 +942,9 @@ def checkCoordinateSystem(self, recover=False, cartesian=False): # Check for differences changed = (IC1 != self.IC) if changed: - logger.info("\x1b[1;94mInternal coordinate system may have changed\x1b[0m") + logger.info("\x1b[1;94mInternal coordinate system may have changed\x1b[0m\n") if self.IC.repr_diff(IC1) != "": - logger.info(self.IC.repr_diff(IC1)) + logger.info(self.IC.repr_diff(IC1)+'\n') # Set current ICs to the new one if changed or recover or cartesian: self.IC = IC1 @@ -999,7 +993,7 @@ def prepareFirstStep(self): # Print initial iteration rms_gradient, max_gradient = self.calcGradNorm() msg = "Step %4i :" % self.Iteration - logger.info(msg + " Gradient = %.3e/%.3e (rms/max) Energy = % .10f" % (rms_gradient, max_gradient, self.E)) + logger.info(msg + " Gradient = %.3e/%.3e (rms/max) Energy = % .10f\n" % (rms_gradient, max_gradient, self.E)) # Initial history self.X_hist = [self.X] self.Gx_hist = [self.gradx] @@ -1029,9 +1023,9 @@ def step(self): v0 = 0.0 if params.verbose: self.IC.Prims.printRotations(self.X) if len(Eig) >= 6: - logger.info("Hessian Eigenvalues: %.5e %.5e %.5e ... %.5e %.5e %.5e" % (Eig[0],Eig[1],Eig[2],Eig[-3],Eig[-2],Eig[-1])) + logger.info("Hessian Eigenvalues: %.5e %.5e %.5e ... %.5e %.5e %.5e\n" % (Eig[0],Eig[1],Eig[2],Eig[-3],Eig[-2],Eig[-1])) else: - logger.info("Hessian Eigenvalues:", ' '.join("%.5e" % i for i in Eig)) + logger.info("Hessian Eigenvalues: " + ' '.join("%.5e" % i for i in Eig) + '\n') # Are we far from constraint satisfaction? self.farConstraints = self.IC.haveConstraints() and self.IC.getConstraintViolation(self.X) > 1e-1 ### OBTAIN AN OPTIMIZATION STEP ### @@ -1042,7 +1036,7 @@ def step(self): inorm = np.linalg.norm(dy) # Cartesian coordinate step size self.cnorm = self.getCartesianNorm(dy) - if params.verbose: logger.info("dy(i): %.4f dy(c) -> target: %.4f -> %.4f" % (inorm, self.cnorm, self.trust)) + if params.verbose: logger.info("dy(i): %.4f dy(c) -> target: %.4f -> %.4f\n" % (inorm, self.cnorm, self.trust)) # If the step is above the trust radius in Cartesian coordinates, then # do the following to reduce the step length: if self.cnorm > 1.1 * self.trust: @@ -1054,31 +1048,31 @@ def step(self): iopt = brent_wiki(froot.evaluate, 0.0, inorm, self.trust, cvg=0.1, obj=froot, verbose=params.verbose) if froot.brentFailed and froot.stored_arg is not None: # If Brent fails but we obtained an IC step that is smaller than the Cartesian trust radius, use it - if params.verbose: logger.info("\x1b[93mUsing stored solution at %.3e\x1b[0m" % froot.stored_val) + if params.verbose: logger.info("\x1b[93mUsing stored solution at %.3e\x1b[0m\n" % froot.stored_val) iopt = froot.stored_arg elif self.IC.bork: # Decrease the target Cartesian step size and try again for i in range(3): froot.target /= 2 - if params.verbose: logger.info("\x1b[93mReducing target to %.3e\x1b[0m" % froot.target) + if params.verbose: logger.info("\x1b[93mReducing target to %.3e\x1b[0m\n" % froot.target) froot.above_flag = True # Stop at any valid step between current target step size and trust radius iopt = brent_wiki(froot.evaluate, 0.0, iopt, froot.target, cvg=0.1, verbose=params.verbose) if not self.IC.bork: break LastForce = self.ForceRebuild self.ForceRebuild = False if self.IC.bork: - logger.info("\x1b[91mInverse iteration for Cartesians failed\x1b[0m") + logger.info("\x1b[91mInverse iteration for Cartesians failed\x1b[0m\n") # This variable is added because IC.bork is unset later. self.ForceRebuild = True else: - if params.verbose: logger.info("\x1b[93mBrent algorithm requires %i evaluations\x1b[0m" % froot.counter) + if params.verbose: logger.info("\x1b[93mBrent algorithm requires %i evaluations\x1b[0m\n" % froot.counter) ##### If IC failed to produce valid Cartesian step, it is "borked" and we need to rebuild it. if self.ForceRebuild: # Force a rebuild of the coordinate system and skip the energy / gradient and evaluation steps. if LastForce: - logger.warning("\x1b[1;91mFailed twice in a row to rebuild the coordinate system; continuing in Cartesian coordinates\x1b[0m") + logger.warning("\x1b[1;91mFailed twice in a row to rebuild the coordinate system; continuing in Cartesian coordinates\x1b[0m\n") self.checkCoordinateSystem(recover=True, cartesian=LastForce) - logger.info("\x1b[1;93mSkipping optimization step\x1b[0m") + logger.info("\x1b[1;93mSkipping optimization step\x1b[0m\n") self.Iteration -= 1 self.state = OPT_STATE.SKIP_EVALUATION return @@ -1140,31 +1134,31 @@ def evaluateStep(self): msg += " Trust = %.3e (%s)" % (self.trust, self.trustprint) msg += " Grad%s = %s%.3e\x1b[0m/%s%.3e\x1b[0m (rms/max)" % ("_T" if self.IC.haveConstraints() else "", "\x1b[92m" if Converged_grms else "\x1b[0m", rms_gradient, "\x1b[92m" if Converged_gmax else "\x1b[0m", max_gradient) # print "Dy.G = %.3f" % Dot, - logger.info(msg + " E (change) = % .10f (%s%+.3e\x1b[0m) Quality = %s%.3f\x1b[0m" % (self.E, "\x1b[91m" if BadStep else ("\x1b[92m" if Converged_energy else "\x1b[0m"), self.E-self.Eprev, "\x1b[91m" if BadStep else "\x1b[0m", Quality)) + logger.info(msg + " E (change) = % .10f (%s%+.3e\x1b[0m) Quality = %s%.3f\x1b[0m" % (self.E, "\x1b[91m" if BadStep else ("\x1b[92m" if Converged_energy else "\x1b[0m"), self.E-self.Eprev, "\x1b[91m" if BadStep else "\x1b[0m", Quality) + "\n") if self.IC is not None and self.IC.haveConstraints(): self.IC.printConstraints(self.X, thre=1e-3) if isinstance(self.IC, PrimitiveInternalCoordinates): - logger.info(self.prim_msg) + logger.info(self.prim_msg + '\n') ### Check convergence criteria ### if Converged_energy and Converged_grms and Converged_drms and Converged_gmax and Converged_dmax and self.conSatisfied: - logger.info("Converged! =D") + logger.info("Converged! =D\n") self.state = OPT_STATE.CONVERGED return if self.Iteration > params.maxiter: - logger.info("Maximum iterations reached (%i); increase --maxiter for more" % params.maxiter) + logger.info("Maximum iterations reached (%i); increase --maxiter for more\n" % params.maxiter) self.state = OPT_STATE.FAILED return if params.qccnv and Converged_grms and (Converged_drms or Converged_energy) and self.conSatisfied: - logger.info("Converged! (Q-Chem style criteria requires grms and either drms or energy)") + logger.info("Converged! (Q-Chem style criteria requires grms and either drms or energy)\n") self.state = OPT_STATE.CONVERGED return if params.molcnv and molpro_converged_gmax and (molpro_converged_dmax or Converged_energy) and self.conSatisfied: - logger.info("Converged! (Molpro style criteria requires gmax and either dmax or energy) This is approximate since convergence checks are done in cartesian coordinates.") + logger.info("Converged! (Molpro style criteria requires gmax and either dmax or energy) This is approximate since convergence checks are done in cartesian coordinates.\n") self.state = OPT_STATE.CONVERGED return @@ -1206,9 +1200,9 @@ def evaluateStep(self): # This is because some systems (e.g. formate) have discontinuities on the # potential surface that can cause an infinite loop if Quality < -1: - if self.trust < self.thre_rj: logger.info("\x1b[93mNot rejecting step - trust below %.3e\x1b[0m" % self.thre_rj) - elif self.E < self.Eprev: logger.info("\x1b[93mNot rejecting step - energy decreases\x1b[0m") - elif self.farConstraints: logger.info("\x1b[93mNot rejecting step - far from constraint satisfaction\x1b[0m") + if self.trust < self.thre_rj: logger.info("\x1b[93mNot rejecting step - trust below %.3e\x1b[0m\n" % self.thre_rj) + elif self.E < self.Eprev: logger.info("\x1b[93mNot rejecting step - energy decreases\x1b[0m\n") + elif self.farConstraints: logger.info("\x1b[93mNot rejecting step - far from constraint satisfaction\x1b[0m\n") # Append steps to history (for rebuilding Hessian) self.X_hist.append(self.X) @@ -1217,16 +1211,16 @@ def evaluateStep(self): ### Rebuild Coordinate System if Necessary ### UpdateHessian = True if self.IC.bork: - logger.info("Failed inverse iteration - checking coordinate system") + logger.info("Failed inverse iteration - checking coordinate system\n") self.checkCoordinateSystem(recover=True) UpdateHessian = False elif self.CoordCounter == (params.check - 1): - logger.info("Checking coordinate system as requested every %i cycles" % params.check) + logger.info("Checking coordinate system as requested every %i cycles\n" % params.check) if self.checkCoordinateSystem(): UpdateHessian = False else: self.CoordCounter += 1 if self.IC.largeRots(): - logger.info("Large rotations - refreshing Rotator reference points and DLC vectors") + logger.info("Large rotations - refreshing Rotator reference points and DLC vectors\n") self.refreshCoordinates() UpdateHessian = False self.G = self.IC.calcGrad(self.X, self.gradx).flatten() @@ -1257,9 +1251,9 @@ def evaluateStep(self): Eig1.sort() if params.verbose: msg += " Eig-ratios: %.5e ... %.5e" % (np.min(Eig1)/np.min(Eig), np.max(Eig1)/np.max(Eig)) - logger.info(msg) + logger.info(msg+'\n') if np.min(Eig1) <= params.epsilon and params.reset: - logger.info("Eigenvalues below %.4e (%.4e) - returning guess" % (params.epsilon, np.min(Eig1))) + logger.info("Eigenvalues below %.4e (%.4e) - returning guess\n" % (params.epsilon, np.min(Eig1))) self.H = self.IC.guess_hessian(self.coords) # Then it's on to the next loop iteration! return @@ -1325,8 +1319,8 @@ def CheckInternalGrad(coords, molecule, IC, engine, dirname, verbose=False): # Initial internal coordinates q0 = IC.calculate(coords) Gq = IC.calcGrad(coords, gradx) - logger.info("-=# Now checking gradient of the energy in internal coordinates vs. finite difference #=-") - logger.info("%20s : %14s %14s %14s" % ('IC Name', 'Analytic', 'Numerical', 'Abs-Diff')) + logger.info("-=# Now checking gradient of the energy in internal coordinates vs. finite difference #=-\n") + logger.info("%20s : %14s %14s %14s\n" % ('IC Name', 'Analytic', 'Numerical', 'Abs-Diff')) h = 1e-3 Gq_f = np.zeros_like(Gq) for i in range(len(q0)): @@ -1338,7 +1332,7 @@ def CheckInternalGrad(coords, molecule, IC, engine, dirname, verbose=False): x1 = IC.newCartesian(coords, dq, verbose) EMinus, _ = engine.calc(x1, dirname) fdiff = (EPlus-EMinus)/(2*h) - logger.info("%20s : % 14.6f % 14.6f % 14.6f" % (IC.Internals[i], Gq[i], fdiff, Gq[i]-fdiff)) + logger.info("%20s : % 14.6f % 14.6f % 14.6f\n" % (IC.Internals[i], Gq[i], fdiff, Gq[i]-fdiff)) Gq_f[i] = fdiff return Gq, Gq_f @@ -1354,9 +1348,9 @@ def CheckInternalHess(coords, molecule, IC, engine, dirname, verbose=False): # Calculate Hessian using finite difference nc = len(coords) Hx = np.zeros((nc, nc), dtype=float) - logger.info("Calculating Cartesian Hessian using finite difference on Cartesian gradients") + logger.info("Calculating Cartesian Hessian using finite difference on Cartesian gradients\n") for i in range(nc): - logger.info(" coordinate %i/%i" % (i+1, nc)) + logger.info(" coordinate %i/%i\n" % (i+1, nc)) coords[i] += h _, gplus = engine.calc(coords, dirname) coords[i] -= 2*h @@ -1372,8 +1366,8 @@ def CheckInternalHess(coords, molecule, IC, engine, dirname, verbose=False): Gq = IC.calcGrad(coords, gradx) Hq_f = np.zeros((len(q0), len(q0)), dtype=float) - logger.info("-=# Now checking Hessian of the energy in internal coordinates using finite difference on gradient #=-") - logger.info("%20s %20s : %14s %14s %14s" % ('IC1 Name', 'IC2 Name', 'Analytic', 'Numerical', 'Abs-Diff')) + logger.info("-=# Now checking Hessian of the energy in internal coordinates using finite difference on gradient #=-\n") + logger.info("%20s %20s : %14s %14s %14s\n" % ('IC1 Name', 'IC2 Name', 'Analytic', 'Numerical', 'Abs-Diff')) h = 1.0e-2 for i in range(len(q0)): dq = np.zeros_like(q0) @@ -1391,7 +1385,7 @@ def CheckInternalHess(coords, molecule, IC, engine, dirname, verbose=False): for j in range(len(q0)): fdiff = fdiffg[j] Hq_f[i, j] = fdiff - logger.info("%20s %20s : % 14.6f % 14.6f % 14.6f" % (IC.Internals[i], IC.Internals[j], Hq[i, j], fdiff, Hq[i, j]-fdiff)) + logger.info("%20s %20s : % 14.6f % 14.6f % 14.6f\n" % (IC.Internals[i], IC.Internals[j], Hq[i, j], fdiff, Hq[i, j]-fdiff)) return Hq, Hq_f @@ -1441,7 +1435,7 @@ def WriteDisplacements(coords, M, IC, dirname, verbose): else: dx = 0.0 x.append((coords+dx).reshape(-1,3) * bohr2ang) - logger.info(i, j, "Displacement (rms/max) = %.5f / %.5f" % (rms_displacement, max_displacement), "(Bork)" if IC.bork else "(Good)") + logger.info("%i %.1f Displacement (rms/max) = %.5f / %.5f %s\n" % (i, j, rms_displacement, max_displacement, "(Bork)" if IC.bork else "(Good)")) M.xyzs = x M.write("%s/ic_%03i.xyz" % (dirname, i)) @@ -1607,7 +1601,6 @@ def run_optimizer(**kwargs): else: logIni = kwargs.get('logIni') logfilename = kwargs.get('prefix') - # Input file for optimization; QC input file or OpenMM .xml file inputf = kwargs.get('input') # Get calculation prefix and temporary directory name @@ -1622,11 +1615,12 @@ def run_optimizer(**kwargs): #| End log file configuration |# #==============================# - logger.info('-=# \x1b[1;94m geomeTRIC started. Version: %s \x1b[0m #=-' % geometric.__version__) - logger.info('geometric-optimize called with the following command line:') - logger.info(' '.join(sys.argv)) + import geometric + logger.info('-=# \x1b[1;94m geomeTRIC started. Version: %s \x1b[0m #=-\n' % geometric.__version__) + logger.info('geometric-optimize called with the following command line:\n') + logger.info(' '.join(sys.argv)+'\n') if backed_up: - logger.info('Backed up existing log file: %s -> %s' % (logfilename, os.path.basename(backed_up))) + logger.info('Backed up existing log file: %s -> %s\n' % (logfilename, os.path.basename(backed_up))) t0 = time.time() params = OptParams(**kwargs) @@ -1644,7 +1638,7 @@ def run_optimizer(**kwargs): # auxiliary functions might (such as writing out displacements) os.makedirs(dirname) else: - logger.info("%s exists ; make sure nothing else is writing to the folder" % dirname) + logger.info("%s exists ; make sure nothing else is writing to the folder\n" % dirname) # TC-specific: Remove existing scratch files in ./run.tmp/scr to avoid confusion for f in ['c0', 'ca0', 'cb0']: if os.path.exists(os.path.join(dirname, 'scr', f)): @@ -1714,10 +1708,11 @@ def run_optimizer(**kwargs): # Print out information about the coordinate system if isinstance(IC, CartesianCoordinates): - logger.info("%i Cartesian coordinates being used" % (3*M.na)) + logger.info("%i Cartesian coordinates being used\n" % (3*M.na)) else: - logger.info("%i internal coordinates being used (instead of %i Cartesians)" % (len(IC.Internals), 3*M.na)) + logger.info("%i internal coordinates being used (instead of %i Cartesians)\n" % (len(IC.Internals), 3*M.na)) logger.info(IC) + logger.info("\n") if Cons is None: # Run a standard geometry optimization @@ -1733,7 +1728,7 @@ def run_optimizer(**kwargs): Mfinal = None for ic, CVal in enumerate(CVals): if len(CVals) > 1: - logger.info("---=== Scan %i/%i : Constrained Optimization ===---" % (ic+1, len(CVals))) + logger.info("---=== Scan %i/%i : Constrained Optimization ===---\n" % (ic+1, len(CVals))) IC = CoordClass(M, build=True, connect=connect, addcart=addcart, constraints=Cons, cvals=CVal, conmethod=params.conmethod) IC.printConstraints(coords, thre=-1) if len(CVals) > 1: @@ -1760,7 +1755,7 @@ def run_optimizer(**kwargs): if len(CVals) > 1: Mfinal.write('scan-final.xyz') print_msg() - logger.info("Time elapsed since start of run_optimizer: %.3f seconds" % (time.time()-t0)) + logger.info("Time elapsed since start of run_optimizer: %.3f seconds\n" % (time.time()-t0)) return progress def main(): diff --git a/geometric/rotate.py b/geometric/rotate.py index e5290bef..d36cf165 100644 --- a/geometric/rotate.py +++ b/geometric/rotate.py @@ -333,7 +333,7 @@ def get_R_der(x, y): RMinus = build_correlation(x, y) x[u, w] += h FDiffR = (RPlus-RMinus)/(2*h) - logger.info(u, w, np.max(np.abs(ADiffR[u, w]-FDiffR))) + logger.info("%i %i %12.6f\n" % (u, w, np.max(np.abs(ADiffR[u, w]-FDiffR)))) return ADiffR def get_F_der(x, y): @@ -398,7 +398,7 @@ def get_F_der(x, y): FMinus = build_F(x, y) x[u, w] += h FDiffF = (FPlus-FMinus)/(2*h) - logger.info(u, w, np.max(np.abs(dF[u, w]-FDiffF))) + logger.info("%i %i %12.6f\n" % (u, w, np.max(np.abs(dF[u, w]-FDiffF)))) return dF def get_q_der(x, y, second=False, fdcheck=False, use_loops=False): @@ -479,7 +479,7 @@ def get_q_der(x, y, second=False, fdcheck=False, use_loops=False): if fdcheck: # If fdcheck = True, then return finite difference derivatives h = 1e-6 - logger.info("-=# Now checking first derivatives of superposition quaternion w/r.t. Cartesians #=-") + logger.info("-=# Now checking first derivatives of superposition quaternion w/r.t. Cartesians #=-\n") FDiffQ = np.zeros((x.shape[0], 3, 4), dtype=float) for u in range(x.shape[0]): for w in range(3): @@ -490,11 +490,11 @@ def get_q_der(x, y, second=False, fdcheck=False, use_loops=False): x[u, w] += h FDiffQ[u, w] = (QPlus-QMinus)/(2*h) maxerr = np.max(np.abs(dq[u, w]-FDiffQ[u, w])) - logger.info("atom %3i %s : maxerr = %.3e %s" % (u, 'xyz'[w], maxerr, 'X' if maxerr > 1e-6 else '')) + logger.info("atom %3i %s : maxerr = %.3e %s\n" % (u, 'xyz'[w], maxerr, 'X' if maxerr > 1e-6 else '')) dq = FDiffQ if second: h = 1.0e-3 - logger.info("-=# Now checking second derivatives of superposition quaternion w/r.t. Cartesians #=-") + logger.info("-=# Now checking second derivatives of superposition quaternion w/r.t. Cartesians #=-\n") Q0 = get_quat(x, y) FDiffQ2 = np.zeros((x.shape[0], 3, y.shape[0], 3, 4), dtype=float) for u in range(x.shape[0]): @@ -523,7 +523,7 @@ def get_q_der(x, y, second=False, fdcheck=False, use_loops=False): FDiffQ2[u, w, a, b] /= (4*h**2) maxerr = np.max(np.abs(dq2[u, w, a, b]-FDiffQ2[u, w, a, b])) if maxerr > 1e-8: - logger.info("atom %3i %s, %3i %s : maxerr = %.3e %s" % (u, 'xyz'[w], a, 'xyz'[b], maxerr, 'X' if maxerr > 1e-6 else '')) + logger.info("atom %3i %s, %3i %s : maxerr = %.3e %s\n" % (u, 'xyz'[w], a, 'xyz'[b], maxerr, 'X' if maxerr > 1e-6 else '')) dq2 = FDiffQ2 if second: return dq, dq2 @@ -633,7 +633,7 @@ def get_expmap_der(x, y, second=False, fdcheck=False, use_loops=False): h = 1e-6 fac, _ = calc_fac_dfac(q[0]) V0 = fac*q[1:] - logger.info("-=# Now checking first derivatives of exponential map w/r.t. quaternion #=-") + logger.info("-=# Now checking first derivatives of exponential map w/r.t. quaternion #=-\n") for p in range(4): # Do backwards difference only, because arccos of q[0] > 1 is undefined q[p] -= h @@ -642,11 +642,11 @@ def get_expmap_der(x, y, second=False, fdcheck=False, use_loops=False): q[p] += h FDiffV = (V0-VMinus)/h maxerr = np.max(np.abs(dvdq[p]-FDiffV)) - logger.info("q %3i : maxerr = %.3e %s" % (i, maxerr, 'X' if maxerr > 1e-6 else '')) + logger.info("q %3i : maxerr = %.3e %s\n" % (i, maxerr, 'X' if maxerr > 1e-6 else '')) # logger.info(i, dvdq[i], FDiffV, np.max(np.abs(dvdq[i]-FDiffV))) if second: h = 1e-5 - logger.info("-=# Now checking second derivatives of exponential map w/r.t. quaternion #=-") + logger.info("-=# Now checking second derivatives of exponential map w/r.t. quaternion #=-\n") def V_(q_): fac_, _ = calc_fac_dfac(q_[0]) return fac_*q_[1:] @@ -673,7 +673,7 @@ def V_(q_): FDiffV2 /= h**2 maxerr = np.max(np.abs(dvdq2[p, r]-FDiffV2)) if maxerr > 1e-7: - logger.info("q %3i %3i : maxerr = %.3e %s" % (p, r, maxerr, 'X' if maxerr > 1e-5 else '')) + logger.info("q %3i %3i : maxerr = %.3e %s\n" % (p, r, maxerr, 'X' if maxerr > 1e-5 else '')) # logger.info("q %3i %3i : analytic %s numerical %s maxerr = %.3e %s" % (i, j, str(dvdq2[i, j]), str(FDiffV2), maxerr, 'X' if maxerr > 1e-4 else '')) # Dimensionality: Number of atoms, number of dimensions (3), number of elements in q (4) @@ -715,7 +715,7 @@ def V_(q_): # print(time.time()-t0) if fdcheck: h = 1e-3 - logger.info("-=# Now checking first derivatives of exponential map w/r.t. Cartesians #=-") + logger.info("-=# Now checking first derivatives of exponential map w/r.t. Cartesians #=-\n") FDiffV = np.zeros((x.shape[0], 3, 3), dtype=float) for u in range(x.shape[0]): for w in range(3): @@ -726,10 +726,10 @@ def V_(q_): x[u, w] += h FDiffV[u, w] = (VPlus-VMinus)/(2*h) maxerr = np.max(np.abs(dvdx[u, w]-FDiffV[u, w])) - logger.info("atom %3i %s : maxerr = %.3e %s" % (u, 'xyz'[w], maxerr, 'X' if maxerr > 1e-6 else '')) + logger.info("atom %3i %s : maxerr = %.3e %s\n" % (u, 'xyz'[w], maxerr, 'X' if maxerr > 1e-6 else '')) dvdx = FDiffV if second: - logger.info("-=# Now checking second derivatives of exponential map w/r.t. Cartesians #=-") + logger.info("-=# Now checking second derivatives of exponential map w/r.t. Cartesians #=-\n") FDiffV2 = np.zeros((x.shape[0], 3, y.shape[0], 3, 3), dtype=float) for u in range(x.shape[0]): for w in range(3): @@ -757,7 +757,7 @@ def V_(q_): FDiffV2[u, w, a, b] /= (4*h**2) maxerr = np.max(np.abs(dvdx2[u, w, a, b]-FDiffV2[u, w, a, b])) if maxerr > 1e-8: - logger.info("atom %3i %s, %3i %s : maxerr = %.3e %s" % (u, 'xyz'[w], a, 'xyz'[b], maxerr, 'X' if maxerr > 1e-6 else '')) + logger.info("atom %3i %s, %3i %s : maxerr = %.3e %s\n" % (u, 'xyz'[w], a, 'xyz'[b], maxerr, 'X' if maxerr > 1e-6 else '')) dvdx2 = FDiffV2 if second: diff --git a/geometric/run_json.py b/geometric/run_json.py index 13404846..37139d6b 100644 --- a/geometric/run_json.py +++ b/geometric/run_json.py @@ -4,6 +4,7 @@ import geometric import json import traceback +import pkg_resources try: from cStringIO import StringIO # Python 2 @@ -11,7 +12,7 @@ from io import StringIO import logging -from .nifty import logger +from .nifty import logger, RawStreamHandler def parse_input_json_dict(in_json_dict): @@ -148,8 +149,13 @@ def make_constraints_string(constraints_dict): def geometric_run_json(in_json_dict): """ Take a input dictionary loaded from json, and return an output dictionary for json """ + # Default logger configuration (prevents extra newline from being printed) + logIni = pkg_resources.resource_filename(geometric.optimize.__name__, 'logJson.ini') + import logging.config + logging.config.fileConfig(logIni,disable_existing_loggers=False) + # Set a temporary logger to capture output - log_stream = logging.StreamHandler(stream=StringIO()) + log_stream = RawStreamHandler(stream=StringIO()) logger.addHandler(log_stream) input_opts = parse_input_json_dict(in_json_dict) @@ -189,10 +195,11 @@ def geometric_run_json(in_json_dict): # Print out information about the coordinate system if isinstance(IC, geometric.internal.CartesianCoordinates): - logger.info("%i Cartesian coordinates being used" % (3 * M.na)) + logger.info("%i Cartesian coordinates being used\n" % (3 * M.na)) else: - logger.info("%i internal coordinates being used (instead of %i Cartesians)" % (len(IC.Internals), 3 * M.na)) + logger.info("%i internal coordinates being used (instead of %i Cartesians)\n" % (len(IC.Internals), 3 * M.na)) logger.info(IC) + logger.info("\n") params = geometric.optimize.OptParams(**input_opts) @@ -208,7 +215,7 @@ def geometric_run_json(in_json_dict): raise RuntimeError("Constraints only work with delocalized internal coordinates") for ic, CVal in enumerate(CVals): if len(CVals) > 1: - logger.info("---=== Scan %i/%i : Constrained Optimization ===---" % (ic + 1, len(CVals))) + logger.info("---=== Scan %i/%i : Constrained Optimization ===---\n" % (ic + 1, len(CVals))) IC = CoordClass(M, build=True, connect=connect, addcart=addcart, constraints=Cons, cvals=CVal) IC.printConstraints(coords, thre=-1) geometric.optimize.Optimize(coords, M, IC, engine, None, params) @@ -238,7 +245,7 @@ def main(): parser.add_argument('in_json', help='Input json file name') parser.add_argument('-o', '--out_json', default='out.json', help='Output Json file name') args = parser.parse_args() - logger.info(' '.join(sys.argv)) + logger.info(' '.join(sys.argv)+"\n") in_json_dict = json.load(open(args.in_json)) out_json_dict = geometric_run_json(in_json_dict) From f18510d7836f5b77c0f1ae16e6594041543cabec Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Sun, 31 Mar 2019 12:23:32 -0700 Subject: [PATCH 2/5] Small changes for Python 2.7 compatibility --- geometric/molecule.py | 2 +- geometric/nifty.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/geometric/molecule.py b/geometric/molecule.py index 2bbc0a9b..19926e53 100644 --- a/geometric/molecule.py +++ b/geometric/molecule.py @@ -215,7 +215,7 @@ def emit(self, record): logger.setLevel(INFO) handler = RawStreamHandler() logger.addHandler(handler) - package = "molecule.py" + package = "LPW-molecule.py" module_name = __name__.replace('.molecule','') diff --git a/geometric/nifty.py b/geometric/nifty.py index f3b2e4e1..ac8b4bc9 100644 --- a/geometric/nifty.py +++ b/geometric/nifty.py @@ -67,6 +67,7 @@ def emit(self, record): message = record.getMessage() self.stream.write(message) self.flush() + class RawFileHandler(FileHandler): """ Exactly like FileHandler, except no newline character is printed at the end of each message. @@ -79,7 +80,9 @@ def __init__(self, *args, **kwargs): def emit(self, record): if self.stream is None: self.stream = self._open() - RawStreamHandler.emit(self, record) + message = record.getMessage() + self.stream.write(message) + self.flush() if "geometric" in __name__: # This ensures logging behavior is consistent with the rest of geomeTRIC @@ -91,7 +94,7 @@ def emit(self, record): logger.setLevel(INFO) handler = RawStreamHandler() logger.addHandler(handler) - package="nifty.py" + package="LPW-nifty.py" try: import bz2 From 7989381c9197a0bfc1baf96a6549f542b90c333f Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Sun, 31 Mar 2019 12:53:45 -0700 Subject: [PATCH 3/5] Small fixes --- geometric/molecule.py | 9 ++++++--- geometric/nifty.py | 5 ++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/geometric/molecule.py b/geometric/molecule.py index 19926e53..83c90bc7 100644 --- a/geometric/molecule.py +++ b/geometric/molecule.py @@ -215,7 +215,10 @@ def emit(self, record): logger.setLevel(INFO) handler = RawStreamHandler() logger.addHandler(handler) - package = "LPW-molecule.py" + if __name__ == "__main__": + package = "LPW-molecule.py" + else: + package = __name__.split('.')[0] module_name = __name__.replace('.molecule','') @@ -1441,7 +1444,7 @@ def __add__(self,other): Sum.Data[key] = self.Data[key] elif diff(self, other, key): for i, j in zip(self.Data[key], other.Data[key]): - logger.info(i, j, i==j) + logger.info("%s %s %s" % (i, j, str(i==j))) logger.error('The data member called %s is not the same for these two objects\n' % key) raise RuntimeError elif key in self.Data: @@ -1483,7 +1486,7 @@ def __iadd__(self,other): if key in ['fnm', 'ftype', 'bonds', 'molecules', 'topology']: pass elif diff(self, other, key): for i, j in zip(self.Data[key], other.Data[key]): - logger.info(i, j, i==j) + logger.info("%s %s %s" % (i, j, str(i==j))) logger.error('The data member called %s is not the same for these two objects\n' % key) raise RuntimeError # Information from the other class is added to this class (if said info doesn't exist.) diff --git a/geometric/nifty.py b/geometric/nifty.py index ac8b4bc9..5bf31f9a 100644 --- a/geometric/nifty.py +++ b/geometric/nifty.py @@ -94,7 +94,10 @@ def emit(self, record): logger.setLevel(INFO) handler = RawStreamHandler() logger.addHandler(handler) - package="LPW-nifty.py" + if __name__ == "__main__": + package = "LPW-nifty.py" + else: + package = __name__.split('.')[0] try: import bz2 From f34c93f4f5690d64c4add9b915efbd294324c6c7 Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Sun, 31 Mar 2019 19:31:16 -0700 Subject: [PATCH 4/5] Simplify AtomContact function --- geometric/molecule.py | 107 ++++++++++++++++++++---------------------- 1 file changed, 51 insertions(+), 56 deletions(-) diff --git a/geometric/molecule.py b/geometric/molecule.py index 83c90bc7..b936abb9 100644 --- a/geometric/molecule.py +++ b/geometric/molecule.py @@ -961,8 +961,9 @@ def AtomContact(xyz, pairs, box=None, displace=False): Parameters ---------- - xyz : np.ndarray - Nx3 array of atom positions + xyz : list or np.ndarray + If a list, must be a list of N_atoms*3 arrays of atom positions. + If an array, must be either a N_atoms*3 (2D) or N_frames*N_atoms*3 (3D). pairs : list List of 2-tuples of atom indices box : np.ndarray, optional @@ -975,45 +976,42 @@ def AtomContact(xyz, pairs, box=None, displace=False): np.ndarray (optional) if displace=True, return a Npairsx3 array of displacement vectors """ + if type(xyz) is list: + xyz = np.array(xyz) # Obtain atom selections for atom pairs parray = np.array(pairs) sel1 = parray[:,0] sel2 = parray[:,1] - xyzpbc = xyz.copy() + if len(xyz.shape) not in (2, 3): + raise RuntimeError("Only provide positions in dimension (traj_length, N_atoms, 3) or (N_atoms, 3)") + single = (len(xyz.shape) == 2) + if single: + xyzpbc = xyz[np.newaxis,:,:].copy() + if box is not None: + box = box[np.newaxis,:].copy() + else: + xyzpbc = xyz.copy() # Minimum image convention: Place all atoms in the box - # [-xbox/2, +xbox/2); [-ybox/2, +ybox/2); [-zbox/2, +zbox/2) + # [0, xbox); [0, ybox); [0, zbox) if box is not None: - xbox = box[0] - ybox = box[1] - zbox = box[2] - while any(xyzpbc[:,0] < -0.5*xbox): - xyzpbc[:,0] += (xyzpbc[:,0] < -0.5*xbox)*xbox - while any(xyzpbc[:,1] < -0.5*ybox): - xyzpbc[:,1] += (xyzpbc[:,1] < -0.5*ybox)*ybox - while any(xyzpbc[:,2] < -0.5*zbox): - xyzpbc[:,2] += (xyzpbc[:,2] < -0.5*zbox)*zbox - while any(xyzpbc[:,0] >= 0.5*xbox): - xyzpbc[:,0] -= (xyzpbc[:,0] >= 0.5*xbox)*xbox - while any(xyzpbc[:,1] >= 0.5*ybox): - xyzpbc[:,1] -= (xyzpbc[:,1] >= 0.5*ybox)*ybox - while any(xyzpbc[:,2] >= 0.5*zbox): - xyzpbc[:,2] -= (xyzpbc[:,2] >= 0.5*zbox)*zbox + box = box[:,np.newaxis,:] + xyzpbc /= box + xyzpbc = xyzpbc % 1.0 # Obtain atom selections for the pairs to be computed # These are typically longer than N but shorter than N^2. - xyzsel1 = xyzpbc[sel1] - xyzsel2 = xyzpbc[sel2] + xyzsel1 = xyzpbc[:,sel1,:] + xyzsel2 = xyzpbc[:,sel2,:] # Calculate xyz displacement dxyz = xyzsel2-xyzsel1 # Apply minimum image convention to displacements if box is not None: - dxyz[:,0] += (dxyz[:,0] < -0.5*xbox)*xbox - dxyz[:,1] += (dxyz[:,1] < -0.5*ybox)*ybox - dxyz[:,2] += (dxyz[:,2] < -0.5*zbox)*zbox - dxyz[:,0] -= (dxyz[:,0] >= 0.5*xbox)*xbox - dxyz[:,1] -= (dxyz[:,1] >= 0.5*ybox)*ybox - dxyz[:,2] -= (dxyz[:,2] >= 0.5*zbox)*zbox - dr2 = np.sum(dxyz**2,axis=1) + dxyz = np.mod(dxyz+0.5, 1.0) - 0.5 + dxyz *= box + dr2 = np.sum(dxyz**2,axis=2) dr = np.sqrt(dr2) + if single: + dr = dr[0] + dxyz = dxyz[0] if displace: return dr, dxyz else: @@ -2158,28 +2156,25 @@ def build_topology(self, force_bonds=True, **kwargs): def distance_matrix(self, pbc=True): """ Obtain distance matrix between all pairs of atoms. """ - AtomIterator = np.ascontiguousarray(np.vstack((np.fromiter(itertools.chain(*[[i]*(self.na-i-1) for i in range(self.na)]),dtype=np.int32), np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) - drij = [] - for sn in range(len(self)): - if hasattr(self, 'boxes') and pbc: - drij.append(AtomContact(self.xyzs[sn],AtomIterator,box=np.array([self.boxes[sn].a, self.boxes[sn].b, self.boxes[sn].c]))) - else: - drij.append(AtomContact(self.xyzs[sn],AtomIterator)) - return AtomIterator, drij + AtomIterator = np.ascontiguousarray(np.vstack((np.fromiter(itertools.chain(*[[i]*(self.na-i-1) for i in range(self.na)]),dtype=np.int32), + np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) + if hasattr(self, 'boxes') and pbc: + boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) + drij = AtomContact(self.xyzs, AtomIterator, box=boxes) + else: + drij = AtomContact(self.xyzs, AtomIterator) + return AtomIterator, list(drij) def distance_displacement(self): """ Obtain distance matrix and displacement vectors between all pairs of atoms. """ - AtomIterator = np.ascontiguousarray(np.vstack((np.fromiter(itertools.chain(*[[i]*(self.na-i-1) for i in range(self.na)]),dtype=np.int32), np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) - drij = [] - dxij = [] - for sn in range(len(self)): - if hasattr(self, 'boxes') and pbc: - drij_i, dxij_i = AtomContact(self.xyzs[sn],AtomIterator,box=np.array([self.boxes[sn].a, self.boxes[sn].b, self.boxes[sn].c]),displace=True) - else: - drij_i, dxij_i = AtomContact(self.xyzs[sn],AtomIterator,box=None,displace=True) - drij.append(drij_i) - dxij.append(dxij_i) - return AtomIterator, drij, dxij + AtomIterator = np.ascontiguousarray(np.vstack((np.fromiter(itertools.chain(*[[i]*(self.na-i-1) for i in range(self.na)]),dtype=np.int32), + np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) + if hasattr(self, 'boxes') and pbc: + boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) + drij, dxij = AtomContact(self.xyzs, AtomIterator, box=boxes, displace=True) + else: + drij, dxij = AtomContact(self.xyzs, AtomIterator, box=None, displace=True) + return AtomIterator, list(drij), list(dxij) def rotate_bond(self, frame, aj, ak, increment=15): """ @@ -2309,21 +2304,21 @@ def find_clashes(self, thre=0.0, pbc=True, groups=None): minDist_frames = [] clashPairs_frames = [] clashDists_frames = [] + if hasattr(self, 'boxes') and pbc: + boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) + drij = AtomContact(self.xyzs, AtomIterator_nb, box=boxes) + else: + drij = AtomContact(self.xyzs, AtomIterator_nb) for frame in range(len(self)): - if hasattr(self, 'boxes') and pbc: - box=np.array([self.boxes[frame].a, self.boxes[frame].b, self.boxes[frame].c]) - drij = AtomContact(self.xyzs[frame],AtomIterator_nb,box=box) - else: - drij = AtomContact(self.xyzs[frame],AtomIterator_nb) - clashPairIdx = np.where(drij < thre)[0] + clashPairIdx = np.where(drij[frame] < thre)[0] clashPairs = AtomIterator_nb[clashPairIdx] - clashDists = drij[clashPairIdx] + clashDists = drij[frame, clashPairIdx] sorter = np.argsort(clashDists) clashPairs = clashPairs[sorter] clashDists = clashDists[sorter] - minIdx = np.argmin(drij) + minIdx = np.argmin(drij[frame]) minPair = AtomIterator_nb[minIdx] - minDist = drij[minIdx] + minDist = drij[frame, minIdx] minPair_frames.append(minPair) minDist_frames.append(minDist) clashPairs_frames.append(clashPairs.copy()) From 4bc53eff15346375f11d776c191bc868dc2f12dc Mon Sep 17 00:00:00 2001 From: Lee-Ping Wang Date: Sun, 31 Mar 2019 19:54:08 -0700 Subject: [PATCH 5/5] Further simplify AtomContact --- geometric/molecule.py | 51 +++++++++++++++++-------------------------- 1 file changed, 20 insertions(+), 31 deletions(-) diff --git a/geometric/molecule.py b/geometric/molecule.py index b936abb9..f8b0c922 100644 --- a/geometric/molecule.py +++ b/geometric/molecule.py @@ -961,41 +961,33 @@ def AtomContact(xyz, pairs, box=None, displace=False): Parameters ---------- - xyz : list or np.ndarray - If a list, must be a list of N_atoms*3 arrays of atom positions. - If an array, must be either a N_atoms*3 (2D) or N_frames*N_atoms*3 (3D). + xyz : np.ndarray + N_frames*N_atoms*3 (3D) array of atomic positions + If you only have a single set of positions, pass in xyz[np.newaxis, :] pairs : list List of 2-tuples of atom indices box : np.ndarray, optional - An array of three numbers (xyz box vectors). + N_frames*3 (2D) array of periodic box vectors + If you only have a single set of positions, pass in box[np.newaxis, :] + displace : bool + If True, also return N_frames*N_pairs*3 array of displacement vectors Returns ------- np.ndarray - A Npairs-length array of minimum image convention distances + N_pairs*N_frames (2D) array of minimum image convention distances np.ndarray (optional) - if displace=True, return a Npairsx3 array of displacement vectors + if displace=True, N_frames*N_pairs*3 array of displacement vectors """ - if type(xyz) is list: - xyz = np.array(xyz) # Obtain atom selections for atom pairs parray = np.array(pairs) sel1 = parray[:,0] sel2 = parray[:,1] - if len(xyz.shape) not in (2, 3): - raise RuntimeError("Only provide positions in dimension (traj_length, N_atoms, 3) or (N_atoms, 3)") - single = (len(xyz.shape) == 2) - if single: - xyzpbc = xyz[np.newaxis,:,:].copy() - if box is not None: - box = box[np.newaxis,:].copy() - else: - xyzpbc = xyz.copy() + xyzpbc = xyz.copy() # Minimum image convention: Place all atoms in the box # [0, xbox); [0, ybox); [0, zbox) if box is not None: - box = box[:,np.newaxis,:] - xyzpbc /= box + xyzpbc /= box[:,np.newaxis,:] xyzpbc = xyzpbc % 1.0 # Obtain atom selections for the pairs to be computed # These are typically longer than N but shorter than N^2. @@ -1006,12 +998,9 @@ def AtomContact(xyz, pairs, box=None, displace=False): # Apply minimum image convention to displacements if box is not None: dxyz = np.mod(dxyz+0.5, 1.0) - 0.5 - dxyz *= box + dxyz *= box[:,np.newaxis,:] dr2 = np.sum(dxyz**2,axis=2) dr = np.sqrt(dr2) - if single: - dr = dr[0] - dxyz = dxyz[0] if displace: return dr, dxyz else: @@ -2072,9 +2061,9 @@ def build_bonds(self): BondThresh = (BT0+BT1) * Fac BondThresh = (BondThresh > mindist) * BondThresh + (BondThresh < mindist) * mindist if hasattr(self, 'boxes') and toppbc: - dxij = AtomContact(self.xyzs[sn], AtomIterator, box=np.array([self.boxes[sn].a, self.boxes[sn].b, self.boxes[sn].c])) + dxij = AtomContact(self.xyzs[sn][np.newaxis, :], AtomIterator, box=np.array([[self.boxes[sn].a, self.boxes[sn].b, self.boxes[sn].c]]))[0] else: - dxij = AtomContact(self.xyzs[sn], AtomIterator) + dxij = AtomContact(self.xyzs[sn][np.newaxis, :], AtomIterator)[0] # Update topology settings with what we learned self.top_settings['toppbc'] = toppbc @@ -2160,9 +2149,9 @@ def distance_matrix(self, pbc=True): np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) if hasattr(self, 'boxes') and pbc: boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) - drij = AtomContact(self.xyzs, AtomIterator, box=boxes) + drij = AtomContact(np.array(self.xyzs), AtomIterator, box=boxes) else: - drij = AtomContact(self.xyzs, AtomIterator) + drij = AtomContact(np.array(self.xyzs), AtomIterator) return AtomIterator, list(drij) def distance_displacement(self): @@ -2171,9 +2160,9 @@ def distance_displacement(self): np.fromiter(itertools.chain(*[range(i+1,self.na) for i in range(self.na)]),dtype=np.int32))).T) if hasattr(self, 'boxes') and pbc: boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) - drij, dxij = AtomContact(self.xyzs, AtomIterator, box=boxes, displace=True) + drij, dxij = AtomContact(np.array(self.xyzs), AtomIterator, box=boxes, displace=True) else: - drij, dxij = AtomContact(self.xyzs, AtomIterator, box=None, displace=True) + drij, dxij = AtomContact(np.array(self.xyzs), AtomIterator, box=None, displace=True) return AtomIterator, list(drij), list(dxij) def rotate_bond(self, frame, aj, ak, increment=15): @@ -2306,9 +2295,9 @@ def find_clashes(self, thre=0.0, pbc=True, groups=None): clashDists_frames = [] if hasattr(self, 'boxes') and pbc: boxes = np.array([[self.boxes[i].a, self.boxes[i].b, self.boxes[i].c] for i in range(len(self))]) - drij = AtomContact(self.xyzs, AtomIterator_nb, box=boxes) + drij = AtomContact(np.array(self.xyzs), AtomIterator_nb, box=boxes) else: - drij = AtomContact(self.xyzs, AtomIterator_nb) + drij = AtomContact(np.array(self.xyzs), AtomIterator_nb) for frame in range(len(self)): clashPairIdx = np.where(drij[frame] < thre)[0] clashPairs = AtomIterator_nb[clashPairIdx]