diff --git a/README.md b/README.md index 526d1e8..ffe2c6a 100644 --- a/README.md +++ b/README.md @@ -103,7 +103,7 @@ inputDict = { "metal_spin": None, # Spin State "full_spin": None, # Assign spin to the full complex (overrides metal_spin) "full_charge": None, # Assign charge to the complex (overrides ligand charges and metal_ox)! - + # Method parameters. "calculator":None, # ASE calculator class input for usage during construction or for optimization. "calculator_kwargs":dict(), # ASE calculator kwargs. diff --git a/architector/io_process_input.py b/architector/io_process_input.py index 27fa439..874d970 100644 --- a/architector/io_process_input.py +++ b/architector/io_process_input.py @@ -467,6 +467,7 @@ def inparse(inputDict): newinpDict = inputDict.copy() u_id = str(uuid.uuid4()) min_n_conformers = 1 + tdebug = newinpDict.get('parameters',{}).get('debug',False) # Adding logging # logging.config.fileConfig('/path/to/logging.conf') if (('core' in inputDict) and (('ligands' in inputDict) or ('ligandList' in inputDict))) or ('mol2string' in inputDict): @@ -584,7 +585,6 @@ def inparse(inputDict): raise ValueError('Need a inputDict["ligands"] list passed') newliglist = [] - tdebug = newinpDict.get('parameters',{}).get('debug',False) if tdebug: print('Starting Ligand Assignment') diff --git a/architector/vibrations_free_energy.py b/architector/vibrations_free_energy.py index 7aeae90..6643b7a 100644 --- a/architector/vibrations_free_energy.py +++ b/architector/vibrations_free_energy.py @@ -96,21 +96,25 @@ def vibration_analysis(atoms, np.dot(centered[i], X[0])*X[j, 1]-np.dot( centered[i], X[1])*X[j, 0]) - ### To implement, qr factorization to normal modes. - + # To implement, qr factorization to normal modes. eig_values, eig_vectors = np.linalg.eigh(mweight_H) - + unit_conversion = ase.units._hbar * ase.units.m / np.sqrt(ase.units._e * ase.units._amu) - - energies = unit_conversion * eig_values.astype(complex)**0.5 # Unit sqrt(eV/Angstroms^2/MW) gives 1/s, convert to eV via frequency - frequencies = energies / ase.units.invcm # Convert to frequency from energy. + # Unit sqrt(eV/Angstroms^2/MW) gives 1/s, convert to eV via frequency + energies = unit_conversion * eig_values.astype(complex)**0.5 + + # Convert to frequency from energy. + frequencies = energies / ase.units.invcm # Modes in columns - mw_normalized = eig_vectors.T # Unitless - md_unnormalized = mw_normalized * mass_weights # Units are 1/sqrt(MW) - norm_factors = 1 / np.linalg.norm(md_unnormalized, axis=0) #units are sqrt(MW) - md_normalized = md_unnormalized * norm_factors # Unitless + mw_normalized = eig_vectors.T # Unitless + # Units are 1/sqrt(MW) + md_unnormalized = mw_normalized * mass_weights + # Units are sqrt(MW) + norm_factors = 1 / np.linalg.norm(md_unnormalized, axis=0) + # Unitless + md_normalized = md_unnormalized * norm_factors # Reshape modes. mw_normalized = mw_normalized.reshape(n_atoms * 3, n_atoms, 3) @@ -119,26 +123,28 @@ def vibration_analysis(atoms, rmasses = norm_factors**2 # units are MW - fconstants = eig_values.astype(complex) * rmasses # units are eV/(Angstroms)^2 to define an R - # kB*T in ase units in units of eV. So plugging into normal mode sampling gives Angstrom. + # Units are eV/(Angstroms)^2 to define an R + fconstants = eig_values.astype(complex) * rmasses + # kB*T in ase units in units of eV. + # So plugging into normal mode sampling gives Angstrom. modes = None - if mode_type == 'direct_eigen_vectors': # Unitless + if mode_type == 'direct_eigen_vectors': # Unitless modes = mw_normalized - elif mode_type == 'mass_weighted_unnormalized': # Mass weighted has units of 1/sqrt(MW) + elif mode_type == 'mass_weighted_unnormalized': # Mass weighted has units of 1/sqrt(MW) modes = md_unnormalized - elif mode_type == 'mass_weighted_normalized': # Unitless + elif mode_type == 'mass_weighted_normalized': # Unitless modes = md_normalized else: raise ValueError('Mode Type requested unknown.') - return energies, modes, fconstants, rmasses, frequencies - + def calc_free_energy(relaxed_atoms, temp=298.15, pressure=101325, - geometry='nonlinear'): + geometry='nonlinear', + ignore_imag_modes=True): """calc_free_energy utility function to calculate free energy of relaxed structures with ASE calculators added. @@ -155,6 +161,8 @@ def calc_free_energy(relaxed_atoms, pressure in pascal, by default 101325 Pa geometry : str, optional 'linear','nonlinear' , by default 'nonlinear' + ignore_imag_modes : bool, optional + default True Returns ------- @@ -174,6 +182,8 @@ def calc_free_energy(relaxed_atoms, potentialenergy=potentialenergy, atoms=relaxed_atoms, geometry=geometry, - symmetrynumber=2, spin=nunpaired/2) + symmetrynumber=2, + spin=nunpaired/2, + ignore_imag_modes=ignore_imag_modes) G = thermo.get_gibbs_energy(temperature=temp, pressure=pressure) - return G, thermo \ No newline at end of file + return G, thermo diff --git a/documentation/tutorials/14-Using_ASE_calcs.ipynb b/documentation/tutorials/14-Using_ASE_calcs.ipynb index d83191b..8c609cb 100644 --- a/documentation/tutorials/14-Using_ASE_calcs.ipynb +++ b/documentation/tutorials/14-Using_ASE_calcs.ipynb @@ -110,8 +110,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## (B) Using arbirtrary ASE calculators to perform energy evaluations.\n", - "\n", + "## (B) How to use example ASE calculator (MACE-MP-0) within Architector complex construction.\n", "Here, we can directly pluging calculators to the full architector search:" ] }, @@ -144,7 +143,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## (C) Using arbirtrary ASE calculators to perform energy evaluations.\n", + "## (C) How to use arbirtrary ASE optimizers within Architector.\n", "\n", "Here, we can directly use other ASE optimizers as well, the default is LBFGSLineSearch. So Let's USE BFGS." ]