diff --git a/README.md b/README.md index 2687b41..b5412b6 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,8 @@ fit_dihedral_with_gomc( atom_type_naming_style='general', gomc_cpu_cores=1, r_squared_min=0.99, - r_squared_rtol=1e-03 + r_squared_rtol=1e-03, + opls_set_k0_zero=True ) import os os.system("cat RB_torsion_k_constants_fit_energy.txt") @@ -53,15 +54,15 @@ The MoSDeF-dihedral-fit documentation can be found [here](https://mosdef-dihedra OPLS-dihedral: -$$OPLS_{Energy} = \frac{f_0}{2}$$ +$$U_{OPLS} = \frac{k_0}{2}$$ -$$+ \frac{f_1}{2} * (1 + cos(\theta)) + \frac{f_2}{2} * (1-cos(2 * \theta))$$ +$$+ \frac{k_1}{2} * (1 + cos(\theta)) + \frac{k_2}{2} * (1-cos(2 * \theta))$$ -$$+ \frac{f_3}{2} * (1 + cos(3 * \theta)) + \frac{f_4}{2} *(1-cos(4 * \theta))$$ +$$+ \frac{k_3}{2} * (1 + cos(3 * \theta)) + \frac{k_4}{2} *(1-cos(4 * \theta))$$ Ryckaert-Bellemans (RB)-torsions: -$$RB_{Energy} = C_0$$ +$$U_{RB} = C_0$$ $$+ C_1 * cos(\psi) + C_2 * cos(\psi)^2$$ @@ -71,7 +72,7 @@ $$\psi = \theta - 180^o$$ Periodic-dihedral: -$$Periodic_{Energy} = K_0 * (1 + cos(n_0*\theta - 90^o))$$ +$$U_{Periodic} = K_0 * (1 + cos(n_0*\theta - 90^o))$$ $$+ K_1 * (1 + cos(n_1*\theta - 180^o)) + K_2 * (1 + cos(n_2*\theta))$$ diff --git a/docs/getting_started/quick_start/fit_a_dihedral_with_Gaussian_and_write_all_files.rst b/docs/getting_started/quick_start/fit_a_dihedral_with_Gaussian_and_write_all_files.rst index 3649d81..786a958 100644 --- a/docs/getting_started/quick_start/fit_a_dihedral_with_Gaussian_and_write_all_files.rst +++ b/docs/getting_started/quick_start/fit_a_dihedral_with_Gaussian_and_write_all_files.rst @@ -83,7 +83,8 @@ Run the dihedral fit to fit to the MM simulations: atom_type_naming_style='general', gomc_cpu_cores=1, r_squared_min=0.99, - r_squared_rtol=1e-03 + r_squared_rtol=1e-03, + opls_force_k0_zero=True ) The most important output files: @@ -195,7 +196,8 @@ Run the dihedral fit to fit to the MM simulations: atom_type_naming_style='general', gomc_cpu_cores=1, r_squared_min=0.99, - r_squared_rtol=5e-03 + r_squared_rtol=5e-03, + opls_force_k0_zero=True ) The most important output files: diff --git a/docs/index.rst b/docs/index.rst index 622275d..d6ae6a0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -41,18 +41,18 @@ the dihedral in **GOMC >= v2.75** and comparing it to the original **Gaussian 16 OPLS-dihedral: .. math:: - OPLS_{Energy} = \frac{f_0}{2} - + \frac{f_1}{2}*(1+cos(\theta)) - + \frac{f_2}{2}*(1-cos(2*\theta)) + U_{OPLS} = \frac{k_0}{2} + + \frac{k_1}{2}*(1+cos(\theta)) + + \frac{k_2}{2}*(1-cos(2*\theta)) .. math:: - + \frac{f_3}{2}*(1+cos(3*\theta)) - + \frac{f_4}{2}*(1-cos(4*\theta)) + + \frac{k_3}{2}*(1+cos(3*\theta)) + + \frac{k_4}{2}*(1-cos(4*\theta)) Ryckaert-Bellemans (RB)-torsions: .. math:: - RB_{Energy} = C_0 + C_1*cos(\psi) + U_{RB} = C_0 + C_1*cos(\psi) + C_2*cos(\psi)^2 + C_3*cos(\psi)^3 + C_4*cos(\psi)^4 @@ -63,7 +63,7 @@ Ryckaert-Bellemans (RB)-torsions: Periodic-dihedral: .. math:: - Periodic_{Energy} = K_0 * (1 + cos(n_0*\theta - 90^o)) + U_{Periodic} = K_0 * (1 + cos(n_0*\theta - 90^o)) .. math:: + K_1 * (1 + cos(n_1*\theta - 180^o)) diff --git a/docs/overview/general_info.rst b/docs/overview/general_info.rst index a6b84ef..95035a7 100644 --- a/docs/overview/general_info.rst +++ b/docs/overview/general_info.rst @@ -44,18 +44,18 @@ the dihedral in **GOMC >= v2.75** and comparing it to the original **Gaussian 16 OPLS-dihedral: .. math:: - OPLS_{Energy} = \frac{f_0}{2} - + \frac{f_1}{2}*(1+cos(\theta)) - + \frac{f_2}{2}*(1-cos(2*\theta)) + U_{OPLS} = \frac{k_0}{2} + + \frac{k_1}{2}*(1+cos(\theta)) + + \frac{k_2}{2}*(1-cos(2*\theta)) .. math:: - + \frac{f_3}{2}*(1+cos(3*\theta)) - + \frac{f_4}{2}*(1-cos(4*\theta)) + + \frac{k_3}{2}*(1+cos(3*\theta)) + + \frac{k_4}{2}*(1-cos(4*\theta)) Ryckaert-Bellemans (RB)-torsions: .. math:: - RB_{Energy} = C_0 + C_1*cos(\psi) + U_{RB} = C_0 + C_1*cos(\psi) + C_2*cos(\psi)^2 + C_3*cos(\psi)^3 + C_4*cos(\psi)^4 @@ -66,7 +66,7 @@ Ryckaert-Bellemans (RB)-torsions: Periodic-dihedral: .. math:: - Periodic_{Energy} = K_0 * (1 + cos(n_0*\theta - 90^o)) + U_{Periodic} = K_0 * (1 + cos(n_0*\theta - 90^o)) .. math:: + K_1 * (1 + cos(n_1*\theta - 180^o)) diff --git a/mosdef_dihedral_fit/dihedral_fit/fit_dihedral_with_gomc.py b/mosdef_dihedral_fit/dihedral_fit/fit_dihedral_with_gomc.py index 4fc2c5a..2bf2643 100755 --- a/mosdef_dihedral_fit/dihedral_fit/fit_dihedral_with_gomc.py +++ b/mosdef_dihedral_fit/dihedral_fit/fit_dihedral_with_gomc.py @@ -38,6 +38,7 @@ def fit_dihedral_with_gomc( gomc_cpu_cores=1, r_squared_min=0.98, r_squared_rtol=2.5e-02, + opls_force_k0_zero=False, ): """Fit the desired dihedral to a MM force field, based on QM data. @@ -67,6 +68,16 @@ def fit_dihedral_with_gomc( and recreated while running this function to ensure only the lasted data is in these folders. + NOTE: The OPLS dihedral equation + + .. math:: + opls_dihedral_n_1 &= 1/2 *( + &= k0 + &= + k1 * (1 + cos(1 * phi)) + &= + k2 * (1 - cos(2 * phi)) + &= + k3 * (1 + cos(3 * phi)) + &= + k4 * (1 - cos(4 * phi)) + &= ) Parameters ---------- @@ -314,6 +325,21 @@ def fit_dihedral_with_gomc( compared to the QM data. NOTE: This value may need adjusted to get the dihedral fit to solve correctly. + opls_force_k0_zero: bool, default=False + The k0 constant is from the equation listed above. + If True, this is force sets the k0 constant in the opls equation to zero (0), + which is the original OPLS form. This means that the dihedral energy fit must + be zero (0) at dihedral angles of (-180 and 180 degrees), which could mean + the dihedral produces both positive (+) and negative energy values (-). + NOTE: Using this option may not allow some dihedrals to be fit properly. + + If False, the k0 constant is allowed to be fitted to a non-zero constant. + This allows the k0 constant to be zero (0) or a scaler value. In this case, + the dihedral's k0 constant is shifted to make the minimum of the dihedral + equation equal to zero (0). This is the more standard overall dihedral form + (i.e. not forcing the equation k0 constant to be zero) as it allows more + dihedrals to be fit properly, and this energy shifting is in almost all newer + dihedral forms now. Outputs ------- @@ -578,8 +604,6 @@ def fit_dihedral_with_gomc( else: raise TypeError(print_error_value) - print("*****************") - print(f"value_j = {str(value_j)}") if isinstance(value_j, list): for int_j in value_j: if not isinstance(int_j, int) or int_j < 0: @@ -716,6 +740,12 @@ def fit_dihedral_with_gomc( f"ERROR: The 'qm_engine' is a {type(qm_engine)}, but it needs to be a str." ) + # check if opls_force_k0_zero' is a bool + if not isinstance(opls_force_k0_zero, bool): + raise TypeError( + "ERROR: Please enter the 'opls_force_k0_zero' as a bool." + ) + # ************************************************************** # ************************************************************** # Use the existing Gaussian (QM) file for their energy data. @@ -764,10 +794,8 @@ def fit_dihedral_with_gomc( # liquid_box_0_length_nm must be a value <= 999.8 nm and an interger value in angstroms <= 9998 Ang liquid_box_0_length_nm = 999.8 - print( - "INFO: Started Building the fragment for the GOMC simulation with dihedral k values = 0" - ) + # Started building the fragment for the GOMC simulation with dihedral k values = 0 box_0_liq = mb.fill_box( compound=[fragment], n_compounds=[1], @@ -779,13 +807,7 @@ def fit_dihedral_with_gomc( seed=seed_no, ) - print( - "INFO: Finished Building the fragment for the GOMC simulation with dihedral k values = 0" - ) - - print( - "INFO: Started Building the Charmm Object for the GOMC simulation with dihedral k values = 0" - ) + # Started building the Charmm object for the GOMC simulation with dihedral k values = 0 # build the charmm object charmm = mf_charmm.Charmm( box_0_liq, @@ -801,22 +823,13 @@ def fit_dihedral_with_gomc( gomc_fix_bonds_angles=None, atom_type_naming_style=atom_type_naming_style, ) - print( - "INFO: Finished Building the Charmm Object for the GOMC simulation with dihedral k values = 0" - ) - # Write the write the FF (.inp), psf, pdb, and GOMC control files - print( - "INFO: Started Writing the PDB, PSF, and FF files for the GOMC simulation with dihedral k values = 0" - ) + # Started writing the PDB, PSF, and FF files for the GOMC simulation with dihedral k values = 0 charmm.write_inp() charmm.write_psf() charmm.write_pdb() - print( - "INFO: Finished Writing the PDB, PSF, and FF files for the GOMC simulation with dihedral k values = 0" - ) # ************************************************************** # make the PDB, PSF and FF files for all the dihedral angles (END) @@ -931,11 +944,7 @@ def fit_dihedral_with_gomc( CrankShaftFreq = 0 SwapFreq = 0 - print("#**********************") - print( - "Started: NVT GOMC control file writing for the GOMC simulation with dihedral k values = 0" - ) - print("#**********************") + # NVT GOMC control file writing for the GOMC simulation with dihedral k values = 0 # calc MC steps for gomc equilb MC_steps = 2 EqSteps = 1 @@ -943,8 +952,6 @@ def fit_dihedral_with_gomc( output_true_list_input = [True, 1] output_false_list_input = [False, 1] - print(f"charmm.combining_rule = {charmm.combining_rule}") - gomc_control.write_gomc_control_file( charmm, f"{gomc_runs_folder_name}/{control_file_name_str}", @@ -995,18 +1002,12 @@ def fit_dihedral_with_gomc( "HistogramFreq": output_false_list_input, "CoordinatesFreq": output_false_list_input, "DCDFreq": output_false_list_input, - "Potential": "VDW", "CBMC_First": 12, "CBMC_Nth": 10, "CBMC_Ang": 50, "CBMC_Dih": 50, }, ) - print("#**********************") - print( - "Completed: NVT GOMC control file written for the GOMC simulation with dihedral k values = 0" - ) - print("#**********************") # ************************************************************** # make the GOMC control file for each dihedral angle (END) @@ -1154,15 +1155,6 @@ def fit_dihedral_with_gomc( for i in GOMC_data_total_energy_kcal_per_mol_list ] - print( - f"GOMC_data_dihedral_degrees_list = {GOMC_data_dihedral_degrees_list}" - ) - print( - f"GOMC_data_total_energy_kcal_per_mol_list = {GOMC_data_total_energy_kcal_per_mol_list}" - ) - print( - f"GOMC_data_total_energy_kcal_per_mol_normalize_list = {GOMC_data_total_energy_kcal_per_mol_normalize_list}" - ) # ********************************* # get GOMC data (END) # ********************************* @@ -1192,17 +1184,6 @@ def fit_dihedral_with_gomc( for i in gaussian_data_total_energy_kcal_per_mol_list ] - print( - f"gaussian_data_dihedral_degrees_list = {gaussian_data_dihedral_degrees_list}" - ) - print( - f"gaussian_data_total_energy_kcal_per_mol_list = {gaussian_data_total_energy_kcal_per_mol_list}" - ) - print( - f"gaussian_data_total_energy_kcal_per_mol_normalize_list = " - f"{gaussian_data_total_energy_kcal_per_mol_normalize_list}" - ) - # get the Gaussian minus GOMC total energy and then it normalized Gaussian_minus_GOMC_data_dihedral_degrees_list = ( GOMC_data_dihedral_degrees_list @@ -1223,13 +1204,6 @@ def fit_dihedral_with_gomc( ) ] - print( - f"Gaussian_minus_GOMC_data_dihedral_degrees_list = {Gaussian_minus_GOMC_data_dihedral_degrees_list}" - ) - print( - f"Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list = \ - {Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list}" - ) # ********************************* # get Gaussian data (END) # ********************************* @@ -1428,16 +1402,6 @@ def fit_dihedral_with_gomc( ) ) - print( - f"sorted_GOMC_data_total_energy_kcal_per_mol_normalize_list = {sorted_GOMC_data_total_energy_kcal_per_mol_normalize_list}" - ) - print( - f"sorted_gaussian_data_total_energy_kcal_per_mol_normalize_list = {sorted_gaussian_data_total_energy_kcal_per_mol_normalize_list}" - ) - print( - f"sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list = {sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list}" - ) - plot_max = int( max( sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list @@ -1646,6 +1610,27 @@ def fit_dihedral_with_gomc( # (END) # ********************************** + # modify 'sorted_const_1_minus_Cos_0_phi_data_lists' based on "opls_force_k0_zero" + # which allows k0=0 or a k0=constant. + # If 'sorted_const_1_minus_Cos_0_phi_data_lists' = all 0s, then k0=0, + # which is changed to 'k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists'. + # If 'sorted_const_1_minus_Cos_0_phi_data_lists' = all 1s, then k0=constant, + # which is changed to 'k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists'. + # It is critical that 'const_1_minus_Cos_0_phi_data' be all (0, 0, .., 0) for k0=0 + # and all (1, 1, .., 1) 'const_1_minus_Cos_0_phi_data' for k0=constant. + k0_forced_to_0_list = [] + k0_is_constant_list = [] + for determine_k0_m in range(len(sorted_const_1_minus_Cos_0_phi_data_lists)): + k0_forced_to_0_list.append(0) + k0_is_constant_list.append(1) + + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists = tuple( + k0_forced_to_0_list + ) + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists = tuple( + k0_is_constant_list + ) + # Run the fitting for only the allowed power types for k_type_i in fit_k_list_allowed: # make the list of k_type_i for fitting in the data @@ -1656,14 +1641,69 @@ def fit_dihedral_with_gomc( k_type_list_i.append(k_type_i) if k_type_i == "1": - parameters, covariance = curve_fit( + # ***** (k=0 -> forced) **** + # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[2] = 0 + parameters_k0_forced_to_0[3] = 0 + parameters_k0_forced_to_0[4] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1675,41 +1715,96 @@ def fit_dihedral_with_gomc( ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[2] = 0 - parameters[3] = 0 - parameters[4] = 0 - - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[2] = 0 + parameters_k0_is_constant[3] = 0 + parameters_k0_is_constant[4] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "2": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[1] = 0 + parameters_k0_forced_to_0[3] = 0 + parameters_k0_forced_to_0[4] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1719,42 +1814,98 @@ def fit_dihedral_with_gomc( sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list ), ) + # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[1] = 0 - parameters[3] = 0 - parameters[4] = 0 - - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[1] = 0 + parameters_k0_is_constant[3] = 0 + parameters_k0_is_constant[4] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "3": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[1] = 0 + parameters_k0_forced_to_0[2] = 0 + parameters_k0_forced_to_0[4] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1765,41 +1916,46 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[1] = 0 - parameters[2] = 0 - parameters[4] = 0 - - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[1] = 0 + parameters_k0_is_constant[2] = 0 + parameters_k0_is_constant[4] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "4": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1810,41 +1966,96 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[1] = 0 - parameters[2] = 0 - parameters[3] = 0 - - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[1] = 0 + parameters_k0_forced_to_0[2] = 0 + parameters_k0_forced_to_0[3] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[1] = 0 + parameters_k0_is_constant[2] = 0 + parameters_k0_is_constant[3] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "1_3": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1855,40 +2066,143 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[2] = 0 - parameters[4] = 0 + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[2] = 0 + parameters_k0_forced_to_0[4] = 0 - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[2] = 0 + parameters_k0_is_constant[4] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "2_4": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[1] = 0 + parameters_k0_forced_to_0[3] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1899,40 +2213,94 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[1] = 0 - parameters[3] = 0 + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[1] = 0 + parameters_k0_is_constant[3] = 0 - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "1_2": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[3] = 0 + parameters_k0_forced_to_0[4] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1943,40 +2311,94 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[3] = 0 - parameters[4] = 0 + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[3] = 0 + parameters_k0_is_constant[4] = 0 - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "3_4": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[1] = 0 + parameters_k0_forced_to_0[2] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -1987,40 +2409,45 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[1] = 0 - parameters[2] = 0 + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[1] = 0 + parameters_k0_is_constant[2] = 0 - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "1_2_3": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -2031,39 +2458,140 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 - parameters[4] = 0 + # and the unused ones are auto-set to 1. (k=0 -> forced) + + parameters_k0_forced_to_0[0] = 0 + parameters_k0_forced_to_0[4] = 0 - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=constant) + parameters_k0_is_constant[4] = 0 + + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) elif k_type_i == "1_2_3_4": - parameters, _ = curve_fit( + # ***** (k=0 -> forced) **** + # # get parameeters and covariance the OPLS dihedral (k=0 -> forced) + parameters_k0_forced_to_0, covariance_k0_forced_to_0 = curve_fit( + mdf_math.opls_dihedral, + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + np.asarray( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list + ), + ) + # fix parameters to zero for the unused values because we want the list length the same, + # and the unused ones are auto-set to 1. (k=0 -> forced) + parameters_k0_forced_to_0[0] = 0 + + # fit the OPLS dihedral (k=0 -> forced) + fit_opls_dihedral_k0_forced_to_0 = mdf_math.opls_dihedral( + ( + np.asarray(k_type_list_i), + np.asarray( + sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list + ), + np.asarray( + k0_forced_to_0_sorted_const_1_minus_Cos_0_phi_data_lists + ), + np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), + np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), + np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), + ), + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + + # ***** (k=constant) **** + # # get parameeters and covariance the OPLS dihedral (k=constant) + parameters_k0_is_constant, covariance_k0_is_constant = curve_fit( mdf_math.opls_dihedral, ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), @@ -2074,27 +2602,28 @@ def fit_dihedral_with_gomc( ), ) # fix parameters to zero for the unused values because we want the list length the same, - # and the unused ones are auto-set to 1. - parameters[0] = 0 + # and the unused ones are auto-set to 1. (k=constant) - # fit the OPLS dihedral - fit_opls_dihedral = mdf_math.opls_dihedral( + # fit the OPLS dihedral (k=constant) + fit_opls_dihedral_k0_is_constant = mdf_math.opls_dihedral( ( np.asarray(k_type_list_i), np.asarray( sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list ), - np.asarray(sorted_const_1_minus_Cos_0_phi_data_lists), + np.asarray( + k0_is_constant_sorted_const_1_minus_Cos_0_phi_data_lists + ), np.asarray(sorted_const_1_plus_Cos_1_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_2_phi_data_lists), np.asarray(sorted_const_1_plus_Cos_3_phi_data_lists), np.asarray(sorted_const_1_minus_Cos_4_phi_data_lists), ), - parameters[0], - parameters[1], - parameters[2], - parameters[3], - parameters[4], + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], ) else: @@ -2102,12 +2631,95 @@ def fit_dihedral_with_gomc( f"ERROR: The {k_type_i} selected in the 'fit_k_list' variable is not a valid selection" ) - # calulate TSS, RSS and R**2 - r_squared = mdf_math.get_r_squared( + # shift dihedral so always at energy = 0 minimum (changing 'fit_opls_dihedral' + C and 'parameters[0]' + C, + # and 'sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list' + C) + theta_spacing_degree = 1 + + # Can not do the min shift for k0_forced_to_0 or you add k0 + fitted_opls_dihedral_small_theta_spacing_k0_forced_to_0 = [] + for theta_i in range(-180, 180, 1): + fitted_opls_dihedral_small_theta_spacing_k0_forced_to_0.append( + mdf_math.opls_dihedral_n_1_2_3_4( + theta_i, + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ) + ) + + # min shift to 0 for k0_is_constant + fitted_opls_dihedral_small_theta_spacing_k0_is_constant = [] + for theta_i in range(-180, 180, 1): + fitted_opls_dihedral_small_theta_spacing_k0_is_constant.append( + mdf_math.opls_dihedral_n_1_2_3_4( + theta_i, + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], + ) + ) + + min_0_fitted_opls_dihedral_small_theta_spacing_k0_is_constant = np.min( + fitted_opls_dihedral_small_theta_spacing_k0_is_constant + ) + + # subtract minimum to k0_is_constant[0] (shift dihedral so always at energy = 0 minimum ) + # remember min times 2 as OPLS_energy = 1/2 (k0+...) + if opls_force_k0_zero is False: + parameters_k0_is_constant[0] -= ( + min_0_fitted_opls_dihedral_small_theta_spacing_k0_is_constant + * 2 + ) + + # subtract minimum to parameters[0] (shift dihedral so always at energy = 0 minimum ) + for dih_energy_m in range(0, len(fit_opls_dihedral_k0_is_constant)): + fit_opls_dihedral_k0_is_constant[ + dih_energy_m + ] -= min_0_fitted_opls_dihedral_small_theta_spacing_k0_is_constant + + # calulate TSS, RSS and R**2 (k0_is_constant) + r_squared_k0_is_constant = mdf_math.get_r_squared( sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list, - fit_opls_dihedral, + fit_opls_dihedral_k0_is_constant, ) + # calulate TSS, RSS and R**2 (k0_forced_to_0) + r_squared_k0_forced_to_0 = mdf_math.get_r_squared( + sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list, + fit_opls_dihedral_k0_forced_to_0, + ) + + # Determine what fit to take 'k0_forced_to_0' or 'k0_is_constant' and + # select the best fit from them + if ( + opls_force_k0_zero is False + and r_squared_k0_is_constant >= r_squared_k0_forced_to_0 + ): + r_squared = r_squared_k0_is_constant + fit_opls_dihedral = fit_opls_dihedral_k0_is_constant + parameters = [ + parameters_k0_is_constant[0], + parameters_k0_is_constant[1], + parameters_k0_is_constant[2], + parameters_k0_is_constant[3], + parameters_k0_is_constant[4], + ] + + else: + r_squared = r_squared_k0_forced_to_0 + fit_opls_dihedral = fit_opls_dihedral_k0_forced_to_0 + parameters = [ + parameters_k0_forced_to_0[0], + parameters_k0_forced_to_0[1], + parameters_k0_forced_to_0[2], + parameters_k0_forced_to_0[3], + parameters_k0_forced_to_0[4], + ] + plt.plot( np.asarray(sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list), np.asarray(fit_opls_dihedral), @@ -2147,7 +2759,7 @@ def fit_dihedral_with_gomc( ncol=2, loc="lower center", fontsize=legend_font_size, - prop={"family": "Arial", "size": legend_font_size}, + prop={"size": legend_font_size}, ) # plt.show() @@ -2469,11 +3081,7 @@ def fit_dihedral_with_gomc( ) output_name_control_fitted_file_name_str = f"output_GOMC_OPLS_fit_{opls_fit_q}_dihedral_coords_{scan_iter_q}.txt" - print("#**********************") - print( - "Started: Writing NVT GOMC control file for the GOMC simulation with fitted dihedral k values." - ) - print("#**********************") + # "Started: Writing NVT GOMC control file for the GOMC simulation with fitted dihedral k values." gomc_control.write_gomc_control_file( charmm, f"{gomc_runs_folder_name}/{control_file_name_fitted_str}", @@ -2524,18 +3132,13 @@ def fit_dihedral_with_gomc( "HistogramFreq": output_false_list_input, "CoordinatesFreq": output_false_list_input, "DCDFreq": output_false_list_input, - "Potential": "VDW", "CBMC_First": 12, "CBMC_Nth": 10, "CBMC_Ang": 50, "CBMC_Dih": 50, }, ) - print("#**********************") - print( - "Completed: NVT GOMC control file written for the GOMC simulation with fitted dihedral k values." - ) - print("#**********************") + # "Completed: NVT GOMC control file written for the GOMC simulation with fitted dihedral k values." # ************************************************************** # make the GOMC control file for each dihedral angle (END) @@ -2545,7 +3148,7 @@ def fit_dihedral_with_gomc( # based on their potential form if charmm.utilized_NB_expression in ["LJ"]: k_constant_units_str = "kcal/mol" - elif charmm.utilized_NB_expression in ["Mie", "Exp6 "]: + elif charmm.utilized_NB_expression in ["Mie", "Exp6"]: k_constant_units_str = "K" else: raise ValueError( @@ -2671,12 +3274,6 @@ def fit_dihedral_with_gomc( GOMC_data_fitted_total_energy_K_list = GOMC_data_fitted_df.loc[ :, "TOTAL" ].tolist() - print( - f"GOMC_data_fitted_dihedral_degrees_list = {GOMC_data_fitted_dihedral_degrees_list}" - ) - print( - f"GOMC_data_fitted_total_energy_K_list = {GOMC_data_fitted_total_energy_K_list}" - ) # convert from Kelvin to kcal/mol normalize so the min value is 0 GOMC_data_fitted_total_energy_kcal_per_mol_list = [ @@ -2688,10 +3285,6 @@ def fit_dihedral_with_gomc( for i in GOMC_data_fitted_total_energy_kcal_per_mol_list ] - print( - f"GOMC_data_fitted_total_energy_kcal_per_mol_normalize_list = {GOMC_data_fitted_total_energy_kcal_per_mol_normalize_list}" - ) - # ********************************* # get GOMC data (END) # ********************************* @@ -2757,11 +3350,6 @@ def fit_dihedral_with_gomc( # Use R**2 (R-squared) to compare the data used to fit the dihedral(s) # via a single or multi-dihedral fit, to the individual fit entered in # GOMC and then recompared to Gaussian and write out the data to a file. - print(f"*****************") - print( - f"GOMC_data_fitted_dihedral_degrees_list === {len(GOMC_data_fitted_dihedral_degrees_list)}" - ) - print(f"*****************") for q_angle in range(0, len(GOMC_data_fitted_dihedral_degrees_list)): if q_angle == 0: # write out the GOMC and Gaussian data in a file @@ -2792,37 +3380,190 @@ def fit_dihedral_with_gomc( f"{str(opls_k_constant_fitted_q_list_kcal_per_mol[4]): <30} " ) - # Compare original fit vs run through GOMC as a validation test case - if opls_fit_data_r_squared_list[ - opls_q - ] >= r_squared_min and not np.isclose( - opls_r_squared_fitted_data_via_gomc_list[opls_q], - opls_fit_data_r_squared_list[opls_q], - rtol=r_squared_rtol, - ): - raise ValueError( - f"ERROR: The calculated R-squared energy values from the fit type " - f"{opls_fit_data_non_zero_k_constants_list[opls_q]} " - f"does not match the validated case for 'r_squared_min' >= " - f"{mdf_math.round_to_sig_figs(r_squared_min,sig_figs=3)}, " - f"within the relative tolerance or 'r_squared_rtol' = " - f"{mdf_math.round_to_sig_figs(r_squared_rtol,sig_figs=3)}. \n" - f"- Fit via the individual or multi-dihedral fit, when " - f"Gaussian minus GOMC with the selected dihedral set to zero \n" - f"--> R-squared = " - f"{mdf_math.round_to_sig_figs(opls_fit_data_r_squared_list[opls_q],sig_figs=3)} \n" - f"- Fit via the validation test case, when " - f"Gaussian minus GOMC with the selected individual dihedral added in GOMC \n" - f"-- >R-squared = " - f"{mdf_math.round_to_sig_figs(opls_r_squared_fitted_data_via_gomc_list[opls_q],sig_figs=3)} \n" - f"The 'r_squared_min' and 'r_squared_rtol' " - f"variables may need to be adjusted, \n" - f"there is likely something wrong with the fitting procedure, the " - f"software parameters need tuned, or there is a bug in the software. \n\n " - f"NOTE: Since the R-squared values are calculated via different parameters, \n" - f"the compared R-squared values could be very different if they are not nearly \n" - f"a perfect fit (R-squared --> ~0.98 to 0.99999999)." + # Compare original fit vs run through EACH GOMC as a validation test case + if opls_q == 0: + opls_fit_acceptable_r_squared_values_list = [] + opls_fit_acceptable_r_squared_values_not_in_rtol_list = [] + opls_fit_data_non_zero_k_constants_not_acceptable_r_squared_values_list = ( + [] + ) + opls_fit_data_r_squared_not_acceptable_r_squared_values_list = [] + opls_r_squared_fitted_data_via_gomc_list_not_acceptable_r_squared_values_list = ( + [] ) + max_opls_fit_Rsquared = None + max_opls_fit_Rsquared_list_index_number = None + + if opls_fit_data_r_squared_list[opls_q] >= r_squared_min: + opls_fit_acceptable_r_squared_values_list.append(1) + + if np.isclose( + opls_r_squared_fitted_data_via_gomc_list[opls_q], + opls_fit_data_r_squared_list[opls_q], + rtol=r_squared_rtol, + ): + opls_fit_acceptable_r_squared_values_not_in_rtol_list.append(1) + + # get maximum value for fit vs run through EACH GOMC as a validation test case + if ( + max_opls_fit_Rsquared is None + and max_opls_fit_Rsquared_list_index_number is None + ): + max_opls_fit_Rsquared = opls_fit_data_r_squared_list[opls_q] + max_opls_fit_Rsquared_list_index_number = opls_q + + elif ( + opls_fit_data_r_squared_list[opls_q] > max_opls_fit_Rsquared + ): + max_opls_fit_Rsquared = opls_fit_data_r_squared_list[opls_q] + max_opls_fit_Rsquared_list_index_number = opls_q + + else: + opls_fit_acceptable_r_squared_values_not_in_rtol_list.append(0) + + opls_fit_data_non_zero_k_constants_not_acceptable_r_squared_values_list.append( + opls_fit_data_non_zero_k_constants_list[opls_q] + ) + opls_fit_data_r_squared_not_acceptable_r_squared_values_list.append( + mdf_math.round_to_sig_figs( + opls_fit_data_r_squared_list[opls_q], sig_figs=8 + ) + ) + opls_r_squared_fitted_data_via_gomc_list_not_acceptable_r_squared_values_list.append( + mdf_math.round_to_sig_figs( + opls_r_squared_fitted_data_via_gomc_list[opls_q], + sig_figs=8, + ) + ) + + # round the 'opls_fit_data_r_squared_list' and 'opls_r_squared_fitted_data_via_gomc_list' for printing + rounded_opls_fit_data_r_squared_list = [ + mdf_math.round_to_sig_figs(r_i, sig_figs=8) + for r_i in opls_fit_data_r_squared_list + ] + rounded_opls_r_squared_fitted_data_via_gomc_list = [ + mdf_math.round_to_sig_figs(r_i, sig_figs=8) + for r_i in opls_r_squared_fitted_data_via_gomc_list + ] + + # combine in pairs for the dihedrals [constants_used, r_squared_fitted, r_squared_rerun_GOMC, tolerance] + rounded_opls_combined_Rsq_tol_list = [] + # Get only ones that meet R^2 objective but not tolerance + meet_r_sq_not_tol_rounded_opls_combined_Rsq_tol_list = [] + for r_n_tol_i in range(0, len(rounded_opls_fit_data_r_squared_list)): + rounded_opls_combined_Rsq_tol_list.append( + [ + opls_fit_data_non_zero_k_constants_list[r_n_tol_i], + rounded_opls_fit_data_r_squared_list[r_n_tol_i], + rounded_opls_r_squared_fitted_data_via_gomc_list[r_n_tol_i], + mdf_math.round_to_sig_figs( + abs( + rounded_opls_fit_data_r_squared_list[r_n_tol_i] + - rounded_opls_r_squared_fitted_data_via_gomc_list[ + r_n_tol_i + ] + ), + sig_figs=8, + ), + ] + ) + + # Get only ones that meet R^2 objective but not tolerance + if opls_fit_data_r_squared_list[r_n_tol_i] >= r_squared_min: + + if not np.isclose( + opls_r_squared_fitted_data_via_gomc_list[r_n_tol_i], + opls_fit_data_r_squared_list[r_n_tol_i], + rtol=r_squared_rtol, + ): + meet_r_sq_not_tol_rounded_opls_combined_Rsq_tol_list.append( + [ + opls_fit_data_non_zero_k_constants_list[r_n_tol_i], + rounded_opls_fit_data_r_squared_list[r_n_tol_i], + rounded_opls_r_squared_fitted_data_via_gomc_list[ + r_n_tol_i + ], + mdf_math.round_to_sig_figs( + abs( + rounded_opls_fit_data_r_squared_list[r_n_tol_i] + - rounded_opls_r_squared_fitted_data_via_gomc_list[ + r_n_tol_i + ] + ), + sig_figs=8, + ), + ] + ) + + # print some Info for the fit, which is the same for all the dihedral angles + print("********************") + print(f"charmm.combining_rule = {charmm.combining_rule}") + print(f"charmm.utilized_NB_expression = {charmm.utilized_NB_expression}") + print("********************") + print( + f"opls_r_squared_fitted_data_via_gomc_list = {opls_r_squared_fitted_data_via_gomc_list}" + ) + print(f"opls_fit_data_r_squared_list = {opls_fit_data_r_squared_list}") + print(f"user set 'r_squared_min' = {r_squared_min}") + print(f"user set 'r_squared_rtol' = {r_squared_rtol}") + print("********************") + print( + f"dihedral_degrees_list = {np.asarray(sorted_Gaussian_minus_GOMC_data_dihedral_degrees_list)}" + ) + print(f"normalized QM - MM total_energy (kcal_per_mol) = ") + print( + f"{np.asarray(sorted_Gaussian_minus_GOMC_data_total_energy_kcal_per_mol_normalized_list)}" + ) + print("********************") + + # Compare original fit vs run through ALL GOMC as a validation test case + if np.sum(opls_fit_acceptable_r_squared_values_not_in_rtol_list) == 0: + raise ValueError( + f"ERROR: The calculated R-squared energy values from the fit types " + f"do not match the validated case for 'r_squared_min' >= " + f"{mdf_math.round_to_sig_figs(r_squared_min,sig_figs=8)}, " + f"within the relative tolerance or 'r_squared_rtol' = " + f"{mdf_math.round_to_sig_figs(r_squared_rtol,sig_figs=8)}. \n" + f"Constants_used = The calculated R-squared energy values from the fit type constants_used." + f"R-squared_fitted = Gaussian minus GOMC with the selected dihedral set to zero \n" + f"R-squared_new_dihedral = Gaussian minus GOMC with the selected individual dihedral added in GOMC \n" + f"Abs(delta) = absolute_value(R-squared_fitted - R-squared_new_dihedral) \n" + f"- The fits for all are shown here for all the dihedral combinations \n" + f"[opls_constants_used, R-squared_fitted, R-squared_new_dihedral, Abs(delta)] \n" + f"{rounded_opls_combined_Rsq_tol_list}, \n" + f"The 'r_squared_min' and 'r_squared_rtol' " + f"variables may need to be adjusted, \n" + f"there is likely something wrong with the fitting procedure, the " + f"software parameters need tuned, or there is a bug in the software. \n\n " + f"NOTE: Since the R-squared values are calculated via different parameters, \n" + f"the compared R-squared values could be very different if they are not nearly \n" + f"a perfect fit (R-squared --> ~0.98 to 0.99999999)." + ) + + elif np.sum( + opls_fit_acceptable_r_squared_values_not_in_rtol_list + ) != np.sum(opls_fit_acceptable_r_squared_values_list): + warn( + f"WARNING: The calculated R-squared energy values from the fit types " + f"do match the validated case for 'r_squared_min' >= " + f"{mdf_math.round_to_sig_figs(r_squared_min,sig_figs=8)}, " + f"but do not fit within the relative tolerance of 'r_squared_rtol' = " + f"{mdf_math.round_to_sig_figs(r_squared_rtol,sig_figs=8)}. \n" + f"Constants_used = The calculated R-squared energy values from the fit type constants_used. \n" + f"R-squared_fitted = Gaussian minus GOMC with the selected dihedral set to zero. \n" + f"R-squared_new_dihedral = Gaussian minus GOMC with the selected individual dihedral added in GOMC. \n" + f"Abs(delta) = Abs(R-squared_fitted - R-squared_new_dihedral). \n" + f"- The fits for all are shown here for all the dihedral combinations \n" + f"[opls_constants_used, R-squared_fitted, R-squared_new_dihedral, Abs(delta)] \n" + f"{meet_r_sq_not_tol_rounded_opls_combined_Rsq_tol_list}. \n" + f"The 'r_squared_min' and 'r_squared_rtol' " + f"variables may need to be adjusted, \n" + f"there is likely something wrong with the fitting procedure, the " + f"software parameters need tuned, or there is a bug in the software. \n\n " + f"NOTE: Since the R-squared values are calculated via different parameters, \n" + f"the compared R-squared values could be very different if they are not nearly \n" + f"a perfect fit (R-squared --> ~0.98 to 0.99999999)." + ) gomc_fitted_gaussian_kcal_per_mol_energy_data_txt_file.close() @@ -2896,7 +3637,7 @@ def fit_dihedral_with_gomc( ncol=2, loc="lower center", fontsize=legend_font_size, - prop={"family": "Arial", "size": legend_font_size}, + prop={"size": legend_font_size}, ) ax2.set_xticks(major_xticks) diff --git a/mosdef_dihedral_fit/tests/files/amber_aa_butane_CT_CT_CT_CT_gmso.xml b/mosdef_dihedral_fit/tests/files/amber_aa_butane_CT_CT_CT_CT_gmso.xml new file mode 100644 index 0000000..5da5b80 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/amber_aa_butane_CT_CT_CT_CT_gmso.xml @@ -0,0 +1,125 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 3 + + + -0.66944 + + + 180 + + + + + + + 3 + + + -0.62760 + + + 180 + + + + + diff --git a/mosdef_dihedral_fit/tests/files/exp6_aa_butane_CT_CT_CT_CT_gmso.xml b/mosdef_dihedral_fit/tests/files/exp6_aa_butane_CT_CT_CT_CT_gmso.xml new file mode 100644 index 0000000..d8c315b --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/exp6_aa_butane_CT_CT_CT_CT_gmso.xml @@ -0,0 +1,141 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/dihedral_molecule_input.xyz b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/dihedral_molecule_input.xyz new file mode 100755 index 0000000..8531bb2 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/dihedral_molecule_input.xyz @@ -0,0 +1,16 @@ +14 +butane_CT_CT_CT_CT +C -2.200 1.212 0.000 +C -1.500 0.000 0.000 +C 0.000 0.000 0.000 +C 0.700 -1.212 -0.000 +H 1.680 -1.025 -0.121 +H -1.813 1.791 0.713 +H -1.967 1.754 -0.813 +H -1.820 -0.485 -0.814 +H -1.820 -0.485 0.814 +H -3.197 1.252 0.096 +H 0.594 -1.558 0.910 +H 0.307 -1.845 -0.496 +H 0.218 0.255 -0.880 +H 0.404 0.678 0.503 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2 b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2 new file mode 100644 index 0000000..a4b499e --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2 @@ -0,0 +1,42 @@ +@MOLECULE +TMP + 14 13 1 0 0 +SMALL +NO_CHARGES +**** +Energy = 0 + +@ATOM + 1 C1 -2.2000 1.2124 0.0000 C 1 TMP 0.000000 + 2 C2 -1.5000 0.0000 0.0000 C 1 TMP 0.000000 + 3 C3 0.0000 0.0000 0.0000 C 1 TMP 0.000000 + 4 C4 0.7000 -1.2124 -0.0000 C 1 TMP 0.000000 + 5 H41 1.6803 -1.0254 -0.1206 H 1 TMP 0.000000 + 6 H11 -1.8131 1.7905 0.7126 H 1 TMP 0.000000 + 7 H12 -1.9666 1.7544 -0.8134 H 1 TMP 0.000000 + 8 H21 -1.8198 -0.4847 -0.8141 H 1 TMP 0.000000 + 9 H22 -1.8198 -0.4847 0.8141 H 1 TMP 0.000000 + 10 H13 -3.1971 1.2516 0.0959 H 1 TMP 0.000000 + 11 H42 0.5936 -1.5578 0.9099 H 1 TMP 0.000000 + 12 H43 0.3070 -1.8446 -0.4962 H 1 TMP 0.000000 + 13 H31 0.2177 0.2554 -0.8800 H 1 TMP 0.000000 + 14 H32 0.4042 0.6781 0.5034 H 1 TMP 0.000000 +@BOND + 1 3 2 1 + 2 3 4 1 + 3 3 13 1 + 4 3 14 1 + 5 2 1 1 + 6 2 8 1 + 7 2 9 1 + 8 1 6 1 + 9 1 7 1 + 10 1 10 1 + 11 4 5 1 + 12 4 11 1 + 13 4 12 1 + +@SUBSTRUCTURE +1 **** 1 TEMP 0 **** **** 0 ROOT + +#generated by VMD diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.pdb b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.pdb new file mode 100644 index 0000000..b5983eb --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.pdb @@ -0,0 +1,16 @@ +CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1 +ATOM 1 C1 TMP X 1 -2.200 1.212 0.000 0.80 0.00 X C +ATOM 2 C2 TMP X 1 -1.500 0.000 0.000 0.80 0.00 X C +ATOM 3 C3 TMP X 1 0.000 0.000 0.000 0.80 0.00 X C +ATOM 4 C4 TMP X 1 0.700 -1.212 -0.000 0.80 0.00 X C +ATOM 5 H41 TMP X 1 1.680 -1.025 -0.121 0.80 0.00 X H +ATOM 6 H11 TMP X 1 -1.813 1.791 0.713 0.80 0.00 X H +ATOM 7 H12 TMP X 1 -1.967 1.754 -0.813 0.80 0.00 X H +ATOM 8 H21 TMP X 1 -1.820 -0.485 -0.814 0.80 0.00 X H +ATOM 9 H22 TMP X 1 -1.820 -0.485 0.814 0.80 0.00 X H +ATOM 10 H13 TMP X 1 -3.197 1.252 0.096 0.80 0.00 X H +ATOM 11 H42 TMP X 1 0.594 -1.558 0.910 0.80 0.00 X H +ATOM 12 H43 TMP X 1 0.307 -1.845 -0.496 0.80 0.00 X H +ATOM 13 H31 TMP X 1 0.218 0.255 -0.880 0.80 0.00 X H +ATOM 14 H32 TMP X 1 0.404 0.678 0.503 0.80 0.00 X H +END diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/.~lock.dihedral.txt# b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/.~lock.dihedral.txt# new file mode 100644 index 0000000..19b6610 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/.~lock.dihedral.txt# @@ -0,0 +1 @@ +,brad,z840-gpu-1,28.07.2024 15:50,file:///home/brad/.config/libreoffice/4; diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral.txt new file mode 100644 index 0000000..bd0ba28 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral.txt @@ -0,0 +1,40 @@ +# Scan of Total Energy +# X-Axis: Scan Coordinate +# Y-Axis: Total Energy (Hartree) +# X Y + -180.0 -157.8299667823 + -170.0 -157.8296078305 + -160.0 -157.8286083701 + -150.0 -157.8271972958 + -140.0 -157.8257264643 + -130.0 -157.8245942859 + -120.0 -157.8241298821 + -110.0 -157.8244742731 + -100.0 -157.8255037773 + -90.0 -157.8268418839 + -80.0 -157.8280112102 + -70.0 -157.8286956505 + -60.0 -157.8287683128 + -50.0 -157.8280936724 + -40.0 -157.8266418577 + -30.0 -157.8246412426 + -20.0 -157.8225262677 + -10.0 -157.8208685196 + 0.0 -157.8202238299 + 10.0 -157.8208681012 + 20.0 -157.8225256601 + 30.0 -157.8246405699 + 40.0 -157.8266416824 + 50.0 -157.8280935628 + 60.0 -157.8287682875 + 70.0 -157.8286957363 + 80.0 -157.8280113682 + 90.0 -157.8268421375 + 100.0 -157.8255039982 + 110.0 -157.8244743633 + 120.0 -157.8241298021 + 130.0 -157.8245940454 + 140.0 -157.8257260663 + 150.0 -157.8271968643 + 160.0 -157.8286080221 + 170.0 -157.8296076205 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_1.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_1.txt new file mode 100644 index 0000000..72d1d1e --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_1.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.264275 1.256039 -0.083201 +2 No Show 2 C -1.500902 -0.059508 0.063771 +3 No Show 3 C 0.013189 0.106939 -0.076323 +4 No Show 4 C 0.776559 -1.208614 0.070608 +5 No Show 5 H 1.846612 -1.057152 -0.033536 +6 No Show 6 H -1.956722 1.974089 0.671989 +7 No Show 7 H -2.087120 1.705324 -1.056471 +8 No Show 8 H -1.854477 -0.767880 -0.682766 +9 No Show 9 H -1.725066 -0.501122 1.032715 +10 No Show 10 H -3.334318 1.104597 0.021072 +11 No Show 11 H 0.599304 -1.657988 1.043819 +12 No Show 12 H 0.469095 -1.926599 -0.684680 +13 No Show 13 H 0.237358 0.548596 -1.045246 +14 No Show 14 H 0.366764 0.815278 0.670247 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_10.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_10.txt new file mode 100644 index 0000000..e69455a --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_10.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.162964 0.936659 -0.476643 +2 No Show 2 C -1.468982 0.109004 0.606437 +3 No Show 3 C 0.064939 0.110903 0.490459 +4 No Show 4 C 0.611472 -1.020699 -0.381587 +5 No Show 5 H 1.693311 -0.969407 -0.459746 +6 No Show 6 H -1.869621 1.980351 -0.410746 +7 No Show 7 H -1.910331 0.585776 -1.473534 +8 No Show 8 H -1.835565 -0.914477 0.582262 +9 No Show 9 H -1.756451 0.509101 1.575503 +10 No Show 10 H -3.243132 0.888808 -0.376339 +11 No Show 11 H 0.354782 -1.988830 0.038552 +12 No Show 12 H 0.207637 -0.981502 -1.389613 +13 No Show 13 H 0.400630 1.070702 0.104997 +14 No Show 14 H 0.500272 0.015609 1.482028 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_11.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_11.txt new file mode 100644 index 0000000..aa53dbd --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_11.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.148216 0.874838 -0.503870 +2 No Show 2 C -1.464022 0.135462 0.647191 +3 No Show 3 C 0.066603 0.098103 0.537568 +4 No Show 4 C 0.591493 -0.969662 -0.423737 +5 No Show 5 H 1.676578 -0.959781 -0.463040 +6 No Show 6 H -1.817135 1.908516 -0.549368 +7 No Show 7 H -1.932437 0.416693 -1.464525 +8 No Show 8 H -1.847532 -0.880338 0.713933 +9 No Show 9 H -1.739281 0.624654 1.578483 +10 No Show 10 H -3.226847 0.878258 -0.378994 +11 No Show 11 H 0.280938 -1.961023 -0.106034 +12 No Show 12 H 0.226975 -0.818135 -1.435409 +13 No Show 13 H 0.432525 1.077645 0.237596 +14 No Show 14 H 0.486353 -0.093233 1.522241 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_12.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_12.txt new file mode 100644 index 0000000..8320392 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_12.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.137775 0.810181 -0.525871 +2 No Show 2 C -1.459564 0.162977 0.682779 +3 No Show 3 C 0.068024 0.082723 0.580064 +4 No Show 4 C 0.576481 -0.914440 -0.462440 +5 No Show 5 H 1.661247 -0.960920 -0.456638 +6 No Show 6 H -1.760460 1.815879 -0.690496 +7 No Show 7 H -1.973749 0.242423 -1.436126 +8 No Show 8 H -1.862706 -0.836147 0.835692 +9 No Show 9 H -1.718246 0.733754 1.571494 +10 No Show 10 H -3.210772 0.880923 -0.375584 +11 No Show 11 H 0.202116 -1.914047 -0.258646 +12 No Show 12 H 0.268045 -0.644401 -1.467310 +13 No Show 13 H 0.466361 1.071942 0.363437 +14 No Show 14 H 0.466988 -0.198850 1.551692 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_13.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_13.txt new file mode 100644 index 0000000..46f2fc0 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_13.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.128277 0.744676 -0.543182 +2 No Show 2 C -1.455467 0.191757 0.714172 +3 No Show 3 C 0.069215 0.064866 0.618936 +4 No Show 4 C 0.563069 -0.857016 -0.497117 +5 No Show 5 H 1.642415 -0.967054 -0.455651 +6 No Show 6 H -1.706620 1.706768 -0.822911 +7 No Show 7 H -2.016283 0.075475 -1.389614 +8 No Show 8 H -1.880533 -0.783094 0.946222 +9 No Show 9 H -1.696393 0.835725 1.556492 +10 No Show 10 H -3.192140 0.886636 -0.379294 +11 No Show 11 H 0.127006 -1.848806 -0.408430 +12 No Show 12 H 0.313138 -0.472210 -1.480111 +13 No Show 13 H 0.501376 1.054345 0.480983 +14 No Show 14 H 0.445480 -0.300066 1.571561 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_14.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_14.txt new file mode 100644 index 0000000..7b361b5 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_14.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.120894 0.683997 -0.556929 +2 No Show 2 C -1.451545 0.224566 0.739307 +3 No Show 3 C 0.069754 0.041278 0.653104 +4 No Show 4 C 0.552381 -0.803164 -0.527287 +5 No Show 5 H 1.624479 -0.964260 -0.468763 +6 No Show 6 H -1.671122 1.599681 -0.932867 +7 No Show 7 H -2.045601 -0.066217 -1.336442 +8 No Show 8 H -1.904103 -0.716536 1.047114 +9 No Show 9 H -1.670053 0.937361 1.530336 +10 No Show 10 H -3.176255 0.879603 -0.393460 +11 No Show 11 H 0.073373 -1.779131 -0.539130 +12 No Show 12 H 0.346709 -0.321783 -1.477078 +13 No Show 13 H 0.540964 1.020838 0.592682 +14 No Show 14 H 0.417895 -0.404223 1.581477 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_15.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_15.txt new file mode 100644 index 0000000..02f6ada --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_15.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.119549 0.629795 -0.567115 +2 No Show 2 C -1.447304 0.261377 0.756597 +3 No Show 3 C 0.068923 0.011516 0.681133 +4 No Show 4 C 0.548332 -0.754550 -0.553034 +5 No Show 5 H 1.614192 -0.950418 -0.486806 +6 No Show 6 H -1.657559 1.502009 -1.022974 +7 No Show 7 H -2.069193 -0.183583 -1.282348 +8 No Show 8 H -1.932006 -0.634754 1.140111 +9 No Show 9 H -1.635607 1.040586 1.490107 +10 No Show 10 H -3.168407 0.860725 -0.407811 +11 No Show 11 H 0.044553 -1.712939 -0.651497 +12 No Show 12 H 0.375303 -0.194137 -1.464995 +13 No Show 13 H 0.584027 0.970342 0.700901 +14 No Show 14 H 0.380274 -0.513954 1.579801 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_16.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_16.txt new file mode 100644 index 0000000..4770698 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_16.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.125580 0.578550 -0.573224 +2 No Show 2 C -1.442333 0.299066 0.766886 +3 No Show 3 C 0.066369 -0.021183 0.702934 +4 No Show 4 C 0.552237 -0.707692 -0.575138 +5 No Show 5 H 1.611687 -0.933374 -0.500367 +6 No Show 6 H -1.657767 1.404847 -1.102124 +7 No Show 7 H -2.100640 -0.288143 -1.224237 +8 No Show 8 H -1.958437 -0.540865 1.227871 +9 No Show 9 H -1.593693 1.142368 1.434648 +10 No Show 10 H -3.167526 0.840328 -0.415816 +11 No Show 11 H 0.030962 -1.644240 -0.755417 +12 No Show 12 H 0.412535 -0.077659 -1.446422 +13 No Show 13 H 0.625149 0.906756 0.808052 +14 No Show 14 H 0.333015 -0.626744 1.564427 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_17.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_17.txt new file mode 100644 index 0000000..3d9a800 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_17.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.138118 0.529859 -0.575302 +2 No Show 2 C -1.436148 0.336551 0.770987 +3 No Show 3 C 0.061704 -0.055527 0.719023 +4 No Show 4 C 0.563226 -0.662200 -0.593627 +5 No Show 5 H 1.616555 -0.912038 -0.509332 +6 No Show 6 H -1.671219 1.309713 -1.170316 +7 No Show 7 H -2.137870 -0.381022 -1.163540 +8 No Show 8 H -1.980092 -0.434447 1.311894 +9 No Show 9 H -1.542936 1.240937 1.362583 +10 No Show 10 H -3.173167 0.817396 -0.417114 +11 No Show 11 H 0.032119 -1.574468 -0.850342 +12 No Show 12 H 0.456128 0.028335 -1.422734 +13 No Show 13 H 0.661287 0.830211 0.916162 +14 No Show 14 H 0.274509 -0.741284 1.533734 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_18.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_18.txt new file mode 100644 index 0000000..756391b --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_18.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.155525 0.478281 -0.572489 +2 No Show 2 C -1.428699 0.370586 0.770752 +3 No Show 3 C 0.055076 -0.087863 0.730168 +4 No Show 4 C 0.579679 -0.612591 -0.609094 +5 No Show 5 H 1.622928 -0.897567 -0.511445 +6 No Show 6 H -1.685126 1.198872 -1.234513 +7 No Show 7 H -2.191772 -0.473929 -1.090439 +8 No Show 8 H -1.991547 -0.320733 1.392496 +9 No Show 9 H -1.487314 1.330434 1.275223 +10 No Show 10 H -3.179450 0.803108 -0.413806 +11 No Show 11 H 0.033952 -1.488530 -0.946312 +12 No Show 12 H 0.518023 0.137722 -1.389811 +13 No Show 13 H 0.687003 0.745911 1.024745 +14 No Show 14 H 0.208747 -0.851685 1.486605 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_19.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_19.txt new file mode 100644 index 0000000..9116368 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_19.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.175672 0.421728 -0.564392 +2 No Show 2 C -1.419902 0.400119 0.767209 +3 No Show 3 C 0.046530 -0.116877 0.737036 +4 No Show 4 C 0.599502 -0.556708 -0.621489 +5 No Show 5 H 1.625512 -0.894510 -0.510357 +6 No Show 6 H -1.694198 1.063132 -1.295243 +7 No Show 7 H -2.263614 -0.568800 -0.998533 +8 No Show 8 H -1.990735 -0.202867 1.467560 +9 No Show 9 H -1.429252 1.406857 1.175036 +10 No Show 10 H -3.181777 0.800581 -0.411482 +11 No Show 11 H 0.030943 -1.377705 -1.045828 +12 No Show 12 H 0.600368 0.254250 -1.342394 +13 No Show 13 H 0.699872 0.656262 1.131160 +14 No Show 14 H 0.138398 -0.953443 1.423803 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_2.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_2.txt new file mode 100644 index 0000000..61aa8f6 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_2.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.266568 1.245079 -0.132260 +2 No Show 2 C -1.496136 -0.047290 0.134837 +3 No Show 3 C 0.019165 0.116863 -0.004967 +4 No Show 4 C 0.771392 -1.213030 0.020759 +5 No Show 5 H 1.842951 -1.061621 -0.066795 +6 No Show 6 H -1.962490 2.033950 0.550544 +7 No Show 7 H -2.092200 1.599948 -1.144027 +8 No Show 8 H -1.837666 -0.812119 -0.560015 +9 No Show 9 H -1.729521 -0.416251 1.131306 +10 No Show 10 H -3.335666 1.098742 -0.012076 +11 No Show 11 H 0.588646 -1.751804 0.946679 +12 No Show 12 H 0.458839 -1.852737 -0.799665 +13 No Show 13 H 0.237733 0.628233 -0.940464 +14 No Show 14 H 0.387521 0.764037 0.788147 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_20.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_20.txt new file mode 100644 index 0000000..6508985 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_20.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.194807 0.365677 -0.552023 +2 No Show 2 C -1.408601 0.428203 0.760280 +3 No Show 3 C 0.034979 -0.145477 0.740641 +4 No Show 4 C 0.618959 -0.499990 -0.629562 +5 No Show 5 H 1.625242 -0.890939 -0.512646 +6 No Show 6 H -1.704984 0.921581 -1.344069 +7 No Show 7 H -2.331607 -0.654424 -0.897682 +8 No Show 8 H -1.979414 -0.080250 1.531569 +9 No Show 9 H -1.366868 1.470038 1.067075 +10 No Show 10 H -3.181764 0.796479 -0.412599 +11 No Show 11 H 0.031238 -1.257787 -1.136138 +12 No Show 12 H 0.680426 0.364754 -1.283197 +13 No Show 13 H 0.700838 0.558985 1.230134 +14 No Show 14 H 0.062335 -1.044832 1.350310 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_21.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_21.txt new file mode 100644 index 0000000..c691273 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_21.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.212541 0.316531 -0.536530 +2 No Show 2 C -1.393889 0.457697 0.748968 +3 No Show 3 C 0.019446 -0.176667 0.741043 +4 No Show 4 C 0.637647 -0.448877 -0.632403 +5 No Show 5 H 1.629150 -0.875932 -0.515899 +6 No Show 6 H -1.726558 0.797929 -1.377841 +7 No Show 7 H -2.380636 -0.723742 -0.800699 +8 No Show 8 H -1.959648 0.046393 1.579661 +9 No Show 9 H -1.297958 1.520998 0.956530 +10 No Show 10 H -3.185765 0.781283 -0.410558 +11 No Show 11 H 0.044813 -1.150628 -1.208444 +12 No Show 12 H 0.741533 0.458982 -1.219946 +13 No Show 13 H 0.691216 0.453257 1.316565 +14 No Show 14 H -0.020837 -1.125206 1.271648 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_22.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_22.txt new file mode 100644 index 0000000..38f3a0e --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_22.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.231825 0.274012 -0.517874 +2 No Show 2 C -1.376562 0.487604 0.732619 +3 No Show 3 C 0.000598 -0.209711 0.737206 +4 No Show 4 C 0.658479 -0.403161 -0.630492 +5 No Show 5 H 1.639179 -0.854576 -0.514686 +6 No Show 6 H -1.759951 0.688416 -1.401730 +7 No Show 7 H -2.419986 -0.779989 -0.705024 +8 No Show 8 H -1.932241 0.171961 1.611008 +9 No Show 9 H -1.225743 1.559385 0.846168 +10 No Show 10 H -3.195020 0.761525 -0.401493 +11 No Show 11 H 0.071845 -1.054227 -1.268903 +12 No Show 12 H 0.793173 0.540581 -1.152563 +13 No Show 13 H 0.671561 0.343671 1.388016 +14 No Show 14 H -0.107537 -1.193469 1.189845 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_23.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_23.txt new file mode 100644 index 0000000..d1f91ad --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_23.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.257361 0.234783 -0.495313 +2 No Show 2 C -1.358834 0.514967 0.710508 +3 No Show 3 C -0.019545 -0.242065 0.727228 +4 No Show 4 C 0.686139 -0.359541 -0.624832 +5 No Show 5 H 1.652940 -0.839350 -0.506991 +6 No Show 6 H -1.804699 0.574578 -1.420145 +7 No Show 7 H -2.469850 -0.826317 -0.599777 +8 No Show 8 H -1.897713 0.289298 1.626665 +9 No Show 9 H -1.157759 1.584523 0.736738 +10 No Show 10 H -3.207158 0.749656 -0.387615 +11 No Show 11 H 0.110801 -0.952305 -1.327181 +12 No Show 12 H 0.856834 0.615378 -1.074714 +13 No Show 13 H 0.642386 0.237348 1.443240 +14 No Show 14 H -0.190216 -1.248923 1.104293 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_24.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_24.txt new file mode 100644 index 0000000..344a8f5 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_24.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.290052 0.199140 -0.468798 +2 No Show 2 C -1.341836 0.539034 0.682155 +3 No Show 3 C -0.039956 -0.273182 0.710267 +4 No Show 4 C 0.721534 -0.318306 -0.615415 +5 No Show 5 H 1.671397 -0.829769 -0.493204 +6 No Show 6 H -1.862929 0.457397 -1.431606 +7 No Show 7 H -2.529679 -0.861237 -0.485577 +8 No Show 8 H -1.858109 0.398340 1.628313 +9 No Show 9 H -1.096527 1.598251 0.626384 +10 No Show 10 H -3.223177 0.745115 -0.369009 +11 No Show 11 H 0.164031 -0.845399 -1.381896 +12 No Show 12 H 0.931923 0.681782 -0.986442 +13 No Show 13 H 0.605952 0.134809 1.483504 +14 No Show 14 H -0.266613 -1.293939 1.013433 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_25.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_25.txt new file mode 100644 index 0000000..ef89072 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_25.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.331033 0.163515 -0.437548 +2 No Show 2 C -1.327869 0.557506 0.647701 +3 No Show 3 C -0.058389 -0.300880 0.685423 +4 No Show 4 C 0.765823 -0.275846 -0.602746 +5 No Show 5 H 1.688679 -0.834437 -0.479748 +6 No Show 6 H -1.931244 0.319233 -1.433911 +7 No Show 7 H -2.610325 -0.883558 -0.352105 +8 No Show 8 H -1.814613 0.496884 1.618087 +9 No Show 9 H -1.048789 1.600965 0.512903 +10 No Show 10 H -3.238406 0.754028 -0.355185 +11 No Show 11 H 0.228098 -0.715954 -1.435798 +12 No Show 12 H 1.030709 0.741523 -0.879258 +13 No Show 13 H 0.563691 0.038773 1.509981 +14 No Show 14 H -0.330374 -1.329711 0.914321 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_26.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_26.txt new file mode 100644 index 0000000..93a7916 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_26.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.376212 0.126723 -0.401644 +2 No Show 2 C -1.317172 0.571149 0.608599 +3 No Show 3 C -0.074381 -0.325456 0.654264 +4 No Show 4 C 0.814918 -0.230964 -0.586659 +5 No Show 5 H 1.699736 -0.850573 -0.476690 +6 No Show 6 H -2.003665 0.156656 -1.420531 +7 No Show 7 H -2.704193 -0.889245 -0.198823 +8 No Show 8 H -1.768605 0.589427 1.597736 +9 No Show 9 H -1.014836 1.594224 0.393961 +10 No Show 10 H -3.249263 0.770592 -0.355531 +11 No Show 11 H 0.297967 -0.558594 -1.482894 +12 No Show 12 H 1.145853 0.791089 -0.750322 +13 No Show 13 H 0.517326 -0.054549 1.525477 +14 No Show 14 H -0.381521 -1.358434 0.805187 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_27.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_27.txt new file mode 100644 index 0000000..da951a9 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_27.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.421100 0.092588 -0.361705 +2 No Show 2 C -1.307990 0.582764 0.565905 +3 No Show 3 C -0.089448 -0.349209 0.618893 +4 No Show 4 C 0.864381 -0.187395 -0.565916 +5 No Show 5 H 1.710511 -0.862556 -0.480641 +6 No Show 6 H -2.079495 -0.004740 -1.387935 +7 No Show 7 H -2.789300 -0.878478 -0.042906 +8 No Show 8 H -1.718555 0.684119 1.567675 +9 No Show 9 H -0.989582 1.579076 0.266973 +10 No Show 10 H -3.260769 0.781068 -0.361287 +11 No Show 11 H 0.374058 -0.396652 -1.512017 +12 No Show 12 H 1.253096 0.825961 -0.612655 +13 No Show 13 H 0.465583 -0.152742 1.533072 +14 No Show 14 H -0.425436 -1.381762 0.684679 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_28.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_28.txt new file mode 100644 index 0000000..477ef7a --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_28.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.467098 0.064880 -0.318206 +2 No Show 2 C -1.300375 0.592323 0.518606 +3 No Show 3 C -0.103687 -0.372427 0.578335 +4 No Show 4 C 0.915610 -0.148902 -0.540040 +5 No Show 5 H 1.729841 -0.864710 -0.478730 +6 No Show 6 H -2.165184 -0.144695 -1.340781 +7 No Show 7 H -2.858969 -0.855617 0.104648 +8 No Show 8 H -1.662598 0.778182 1.526600 +9 No Show 9 H -0.973907 1.555398 0.133399 +10 No Show 10 H -3.279655 0.784141 -0.357301 +11 No Show 11 H 0.462519 -0.250975 -1.522396 +12 No Show 12 H 1.344122 0.847130 -0.476952 +13 No Show 13 H 0.406371 -0.253523 1.530965 +14 No Show 14 H -0.461037 -1.399163 0.553989 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_29.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_29.txt new file mode 100644 index 0000000..9554779 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_29.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.514553 0.045885 -0.271583 +2 No Show 2 C -1.294457 0.599176 0.466343 +3 No Show 3 C -0.117038 -0.394602 0.532036 +4 No Show 4 C 0.968935 -0.117803 -0.508803 +5 No Show 5 H 1.765135 -0.854167 -0.453847 +6 No Show 6 H -2.265576 -0.245262 -1.288782 +7 No Show 7 H -2.907056 -0.831065 0.234880 +8 No Show 8 H -1.599601 0.869594 1.473754 +9 No Show 9 H -0.968570 1.522638 -0.005763 +10 No Show 10 H -3.310675 0.782409 -0.325507 +11 No Show 11 H 0.567127 -0.141711 -1.518471 +12 No Show 12 H 1.411307 0.861950 -0.353898 +13 No Show 13 H 0.338412 -0.355170 1.518000 +14 No Show 14 H -0.487440 -1.409830 0.413781 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_3.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_3.txt new file mode 100644 index 0000000..fa9b712 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_3.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.265861 1.227670 -0.180649 +2 No Show 2 C -1.492584 -0.032943 0.204531 +3 No Show 3 C 0.026186 0.124312 0.065536 +4 No Show 4 C 0.763209 -1.211032 -0.029974 +5 No Show 5 H 1.836457 -1.064714 -0.105097 +6 No Show 6 H -1.966061 2.077672 0.426790 +7 No Show 7 H -2.088815 1.487306 -1.220328 +8 No Show 8 H -1.825488 -0.850721 -0.431720 +9 No Show 9 H -1.737754 -0.324109 1.223298 +10 No Show 10 H -3.334924 1.089977 -0.050045 +11 No Show 11 H 0.575045 -1.830917 0.842781 +12 No Show 12 H 0.441751 -1.768343 -0.905050 +13 No Show 13 H 0.243601 0.704032 -0.829481 +14 No Show 14 H 0.411235 0.703808 0.901415 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_30.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_30.txt new file mode 100644 index 0000000..27e2e97 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_30.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.560099 0.032261 -0.222055 +2 No Show 2 C -1.290058 0.602865 0.410535 +3 No Show 3 C -0.129477 -0.414866 0.481220 +4 No Show 4 C 1.020921 -0.090903 -0.472927 +5 No Show 5 H 1.806345 -0.838387 -0.412878 +6 No Show 6 H -2.371583 -0.323040 -1.231823 +7 No Show 7 H -2.939227 -0.805574 0.356021 +8 No Show 8 H -1.530702 0.956856 1.409171 +9 No Show 9 H -0.974125 1.480265 -0.148090 +10 No Show 10 H -3.345071 0.780676 -0.276115 +11 No Show 11 H 0.678535 -0.052799 -1.503814 +12 No Show 12 H 1.461442 0.873493 -0.236793 +13 No Show 13 H 0.262802 -0.456267 1.493592 +14 No Show 14 H -0.503752 -1.412536 0.266097 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_31.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_31.txt new file mode 100644 index 0000000..eed1d98 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_31.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.600506 0.022115 -0.170302 +2 No Show 2 C -1.286349 0.603671 0.352083 +3 No Show 3 C -0.141686 -0.433194 0.426942 +4 No Show 4 C 1.068229 -0.066530 -0.433153 +5 No Show 5 H 1.845937 -0.821166 -0.362961 +6 No Show 6 H -2.475784 -0.388716 -1.168600 +7 No Show 7 H -2.957377 -0.777831 0.472418 +8 No Show 8 H -1.457653 1.038245 1.332905 +9 No Show 9 H -0.988928 1.429013 -0.290060 +10 No Show 10 H -3.376525 0.780233 -0.218097 +11 No Show 11 H 0.789349 0.026509 -1.479339 +12 No Show 12 H 1.497068 0.881775 -0.121424 +13 No Show 13 H 0.181244 -0.555197 1.457012 +14 No Show 14 H -0.511069 -1.406882 0.114717 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_32.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_32.txt new file mode 100644 index 0000000..25c1a71 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_32.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.634099 0.014713 -0.117261 +2 No Show 2 C -1.282884 0.601894 0.291144 +3 No Show 3 C -0.154077 -0.449817 0.369507 +4 No Show 4 C 1.109050 -0.044228 -0.390338 +5 No Show 5 H 1.880053 -0.804568 -0.308894 +6 No Show 6 H -2.574534 -0.447941 -1.098462 +7 No Show 7 H -2.962139 -0.746093 0.585911 +8 No Show 8 H -1.383176 1.111163 1.245608 +9 No Show 9 H -1.011297 1.370948 -0.428177 +10 No Show 10 H -3.401969 0.781511 -0.157160 +11 No Show 11 H 0.895932 0.101879 -1.445604 +12 No Show 12 H 1.519118 0.885678 -0.005624 +13 No Show 13 H 0.096464 -0.649350 1.407694 +14 No Show 14 H -0.510492 -1.393744 -0.036197 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_33.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_33.txt new file mode 100644 index 0000000..5f9b2bf --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_33.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.660350 0.009815 -0.063890 +2 No Show 2 C -1.279804 0.597847 0.227557 +3 No Show 3 C -0.166527 -0.465084 0.308837 +4 No Show 4 C 1.142716 -0.024047 -0.345382 +5 No Show 5 H 1.907036 -0.789989 -0.254447 +6 No Show 6 H -2.666174 -0.504177 -1.020789 +7 No Show 7 H -2.953995 -0.708295 0.697284 +8 No Show 8 H -1.310617 1.173547 1.148782 +9 No Show 9 H -1.040086 1.308961 -0.559993 +10 No Show 10 H -3.420350 0.784664 -0.097522 +11 No Show 11 H 0.996656 0.176835 -1.402826 +12 No Show 12 H 1.528233 0.883461 0.111896 +13 No Show 13 H 0.011946 -0.736382 1.345940 +14 No Show 14 H -0.502733 -1.375113 -0.183296 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_34.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_34.txt new file mode 100644 index 0000000..bd01502 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_34.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.679404 0.007369 -0.010952 +2 No Show 2 C -1.277704 0.591818 0.161224 +3 No Show 3 C -0.178455 -0.479314 0.244839 +4 No Show 4 C 1.169255 -0.006169 -0.299037 +5 No Show 5 H 1.926509 -0.778464 -0.203402 +6 No Show 6 H -2.749463 -0.560108 -0.934369 +7 No Show 7 H -2.933745 -0.661857 0.806924 +8 No Show 8 H -1.243143 1.224788 1.044423 +9 No Show 9 H -1.074848 1.245875 -0.684405 +10 No Show 10 H -3.431871 0.789531 -0.043112 +11 No Show 11 H 1.090409 0.254330 -1.350494 +12 No Show 12 H 1.525317 0.872760 0.232165 +13 No Show 13 H -0.068900 -0.815181 1.272981 +14 No Show 14 H -0.488008 -1.353336 -0.324627 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_35.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_35.txt new file mode 100644 index 0000000..86a07a6 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_35.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.691652 0.007493 0.041109 +2 No Show 2 C -1.277289 0.584045 0.092401 +3 No Show 3 C -0.189119 -0.492671 0.177722 +4 No Show 4 C 1.188994 0.009149 -0.251748 +5 No Show 5 H 1.939221 -0.770207 -0.158603 +6 No Show 6 H -2.823038 -0.616816 -0.837906 +7 No Show 7 H -2.902072 -0.604581 0.914294 +8 No Show 8 H -1.182908 1.266003 0.934304 +9 No Show 9 H -1.115626 1.183493 -0.801427 +10 No Show 10 H -3.437696 0.795473 0.003438 +11 No Show 11 H 1.176007 0.335818 -1.287513 +12 No Show 12 H 1.511019 0.851308 0.355240 +13 No Show 13 H -0.143657 -0.886328 1.190512 +14 No Show 14 H -0.466235 -1.330136 -0.459662 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_36.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_36.txt new file mode 100644 index 0000000..ca0f454 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_36.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.697418 0.010448 0.092123 +2 No Show 2 C -1.279065 0.574764 0.021786 +3 No Show 3 C -0.197914 -0.505182 0.108142 +4 No Show 4 C 1.202234 0.021602 -0.203640 +5 No Show 5 H 1.946676 -0.764378 -0.120756 +6 No Show 6 H -2.885427 -0.673127 -0.730912 +7 No Show 7 H -2.859245 -0.535771 1.017574 +8 No Show 8 H -1.130805 1.299515 0.819474 +9 No Show 9 H -1.162629 1.122219 -0.911573 +10 No Show 10 H -3.439395 0.801511 0.041936 +11 No Show 11 H 1.252094 0.420360 -1.212855 +12 No Show 12 H 1.485352 0.817883 0.479576 +13 No Show 13 H -0.211238 -0.951742 1.100089 +14 No Show 14 H -0.437272 -1.306061 -0.588798 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_4.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_4.txt new file mode 100644 index 0000000..860f836 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_4.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.261857 1.203954 -0.228421 +2 No Show 2 C -1.490046 -0.016726 0.271818 +3 No Show 3 C 0.033902 0.129229 0.134190 +4 No Show 4 C 0.751712 -1.202761 -0.081549 +5 No Show 5 H 1.826506 -1.065055 -0.151355 +6 No Show 6 H -1.968284 2.104917 0.304074 +7 No Show 7 H -2.075626 1.370972 -1.285332 +8 No Show 8 H -1.818141 -0.884199 -0.297280 +9 No Show 9 H -1.748275 -0.222894 1.307529 +10 No Show 10 H -3.331867 1.076109 -0.095182 +11 No Show 11 H 0.559890 -1.893982 0.735085 +12 No Show 12 H 0.416538 -1.676796 -0.999473 +13 No Show 13 H 0.255265 0.776694 -0.711869 +14 No Show 14 H 0.436282 0.632536 1.009778 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_5.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_5.txt new file mode 100644 index 0000000..18eeb1e --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_5.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.254126 1.173976 -0.275466 +2 No Show 2 C -1.488060 0.001007 0.336026 +3 No Show 3 C 0.041745 0.131759 0.200320 +4 No Show 4 C 0.736493 -1.188217 -0.133783 +5 No Show 5 H 1.812504 -1.060895 -0.205268 +6 No Show 6 H -1.968281 2.116625 0.183883 +7 No Show 7 H -2.053187 1.252503 -1.340047 +8 No Show 8 H -1.815159 -0.912221 -0.156317 +9 No Show 9 H -1.759089 -0.111781 1.382380 +10 No Show 10 H -3.325822 1.055548 -0.146765 +11 No Show 11 H 0.542550 -1.941427 0.625400 +12 No Show 12 H 0.383635 -1.579897 -1.083495 +13 No Show 13 H 0.272323 0.846028 -0.587104 +14 No Show 14 H 0.460470 0.548991 1.112252 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_6.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_6.txt new file mode 100644 index 0000000..138675e --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_6.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.242241 1.137718 -0.321354 +2 No Show 2 C -1.485950 0.019919 0.396914 +3 No Show 3 C 0.049008 0.132166 0.263696 +4 No Show 4 C 0.717193 -1.167245 -0.186215 +5 No Show 5 H 1.794051 -1.051044 -0.264128 +6 No Show 6 H -1.964101 2.114151 0.066076 +7 No Show 7 H -2.023414 1.131490 -1.385510 +8 No Show 8 H -1.815410 -0.933499 -0.009383 +9 No Show 9 H -1.768012 0.008184 1.446065 +10 No Show 10 H -3.315969 1.027967 -0.201942 +11 No Show 11 H 0.521104 -1.974542 0.514290 +12 No Show 12 H 0.344777 -1.477611 -1.158496 +13 No Show 13 H 0.293609 0.910663 -0.455149 +14 No Show 14 H 0.481353 0.453682 1.207157 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_7.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_7.txt new file mode 100644 index 0000000..caf7c49 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_7.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.226041 1.095291 -0.365340 +2 No Show 2 C -1.483064 0.039889 0.454545 +3 No Show 3 C 0.055049 0.130598 0.324439 +4 No Show 4 C 0.693763 -1.139720 -0.238086 +5 No Show 5 H 1.771227 -1.035273 -0.324088 +6 No Show 6 H -1.953658 2.098769 -0.050380 +7 No Show 7 H -1.988980 1.006388 -1.422144 +8 No Show 8 H -1.817669 -0.946054 0.141569 +9 No Show 9 H -1.773129 0.134035 1.497239 +10 No Show 10 H -3.301820 0.994313 -0.257010 +11 No Show 11 H 0.493368 -1.994804 0.401471 +12 No Show 12 H 0.302511 -1.368656 -1.225733 +13 No Show 13 H 0.317667 0.968152 -0.317070 +14 No Show 14 H 0.496773 0.349071 1.292612 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_8.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_8.txt new file mode 100644 index 0000000..0e4ad12 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_8.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.205995 1.047125 -0.406495 +2 No Show 2 C -1.479082 0.061105 0.508993 +3 No Show 3 C 0.059568 0.126907 0.382719 +4 No Show 4 C 0.666815 -1.105784 -0.288460 +5 No Show 5 H 1.744932 -1.014415 -0.380899 +6 No Show 6 H -1.935377 2.071442 -0.166870 +7 No Show 7 H -1.953430 0.875404 -1.449630 +8 No Show 8 H -1.821346 -0.947849 0.293175 +9 No Show 9 H -1.773216 0.261733 1.535505 +10 No Show 10 H -3.283662 0.956683 -0.308083 +11 No Show 11 H 0.457605 -2.003515 0.286162 +12 No Show 12 H 0.260347 -1.251307 -1.285965 +13 No Show 13 H 0.343468 1.015552 -0.175331 +14 No Show 14 H 0.505369 0.238915 1.367204 diff --git a/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_9.txt b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_9.txt new file mode 100644 index 0000000..0a947f3 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/gaussian_style_output_files/CT_CT_CT_CT/output/dihedral_coords_position_9.txt @@ -0,0 +1,15 @@ +Row Highlight Display Tag Symbol X Y Z +1 No Show 1 C -2.183762 0.993989 -0.443891 +2 No Show 2 C -1.474193 0.084003 0.559965 +3 No Show 3 C 0.062719 0.120583 0.438368 +4 No Show 4 C 0.638142 -1.065925 -0.336473 +5 No Show 5 H 1.717476 -0.990756 -0.429062 +6 No Show 6 H -1.907988 2.032488 -0.285539 +7 No Show 7 H -1.923041 0.736471 -1.467200 +8 No Show 8 H -1.826898 -0.937465 0.441322 +9 No Show 9 H -1.767765 0.387540 1.561363 +10 No Show 10 H -3.263024 0.918969 -0.350335 +11 No Show 11 H 0.412250 -2.001593 0.166622 +12 No Show 12 H 0.224561 -1.123504 -1.340011 +13 No Show 13 H 0.370894 1.050256 -0.033448 +14 No Show 14 H 0.506627 0.126942 1.430346 diff --git a/mosdef_dihedral_fit/tests/files/mie_ua_butane_CT_CT_CT_CT_gmso.xml b/mosdef_dihedral_fit/tests/files/mie_ua_butane_CT_CT_CT_CT_gmso.xml new file mode 100644 index 0000000..9b2939c --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/mie_ua_butane_CT_CT_CT_CT_gmso.xml @@ -0,0 +1,227 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + + + 0 + 1 + 2 + 3 + + + 0 + 0 + 0 + 0 + + + 90 + 180 + 0.0 + 180 + + + + + diff --git a/mosdef_dihedral_fit/tests/files/trappe_ua_butane_CT_CT_CT_CT_gmso.xml b/mosdef_dihedral_fit/tests/files/trappe_ua_butane_CT_CT_CT_CT_gmso.xml new file mode 100644 index 0000000..2543752 --- /dev/null +++ b/mosdef_dihedral_fit/tests/files/trappe_ua_butane_CT_CT_CT_CT_gmso.xml @@ -0,0 +1,160 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/mosdef_dihedral_fit/tests/test_file_read_write.py b/mosdef_dihedral_fit/tests/test_file_read_write.py index 325ad39..13478d1 100644 --- a/mosdef_dihedral_fit/tests/test_file_read_write.py +++ b/mosdef_dihedral_fit/tests/test_file_read_write.py @@ -329,7 +329,7 @@ def test_get_matching_dihedral_info_and_opls_fitting_data(self): assert np.shape(degreesList) == (37, 9) opls_paramsList = out[4] expected_opls_params = [ - 0, + 1, 8.999999999999845, 8.99999999999938, 2.780442542871242e-12, diff --git a/mosdef_dihedral_fit/tests/test_fit_dihedral_with_gomc.py b/mosdef_dihedral_fit/tests/test_fit_dihedral_with_gomc.py index 1ff7394..51f4f5b 100644 --- a/mosdef_dihedral_fit/tests/test_fit_dihedral_with_gomc.py +++ b/mosdef_dihedral_fit/tests/test_fit_dihedral_with_gomc.py @@ -17,7 +17,7 @@ # User changable variable, as it needs to be run locally # Examples -# gomc_binary_directory = "/Users/brad/Programs/GOMC/GOMC_2_75/bin" +# gomc_binary_directory = "/home/brad/Programs/GOMC/GOMC-2.75a/bin" # gomc_binary_directory = "/Users/calcraven/Documents/Vanderbilt/Research/MoSDeF/Dihedral_Fitter/GOMC/bin" @@ -50,6 +50,7 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) assert ( @@ -191,7 +192,7 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC(self): "k4_kcal_per_mol", "r_squared", ], - [str(3), 0, 0, 0, 0.314003748427, 0, 0.998791145574], + [str(3), 0, 0, 0, 0.31400374842, 0, 0.998791145574], ] if i == 0: @@ -245,10 +246,10 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC(self): ], [ str("0_3"), - 0.314003748427, + 0.31400374842, 0, 0, - -0.157001874214, + -0.15700187421, 0, 0, 0, @@ -306,10 +307,10 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC(self): ], [ str("0_3"), - 0.157001874212, - 0.471005622638, + 0.15700187421, + 0.47100562263, 0, - -0.62800749685, + -0.62800749684, 0, 0, 0.998791145574, @@ -362,6 +363,7 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC_with_2_log_files( gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) assert ( @@ -503,7 +505,7 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC_with_2_log_files( "k4_kcal_per_mol", "r_squared", ], - [str(3), 0, 0, 0, 0.314003748425, 0, 0.998791145574], + [str(3), 0, 0, 0, 0.31400374842, 0, 0.998791145574], ] if i == 0: @@ -557,10 +559,10 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC_with_2_log_files( ], [ str("0_3"), - 0.314003748425, + 0.31400374842, 0, 0, - -0.157001874212, + -0.15700187421, 0, 0, 0, @@ -618,10 +620,10 @@ def test_gaussian_log_file_fit_oplsaa_fit_ethane_HC_CT_CT_HC_with_2_log_files( ], [ str("0_3"), - 0.157001874212, - 0.471005622638, + 0.15700187421, + 0.47100562263, 0, - -0.62800749685, + -0.62800749684, 0, 0, 0.998791145574, @@ -671,6 +673,7 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) assert ( @@ -809,16 +812,16 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( "k4_kcal_per_mol", "r_squared", ], - [str("1"), 0, 2.0040856888, 0, 0, 0, -0.371866278201], - [str("2"), 0, 0, 2.13847194795, 0, 0, 0.264787158812], - [str("3"), 0, 0, 0, 1.79252468631, 0, -1.29046011346], + [str("1"), 0, 2.00408568934, 0, 0, 0, -0.371866278201], + [str("2"), 0, 0, 2.13847194875, 0, 0, 0.264787158812], + [str("3"), 0, 0, 0, 1.79252468736, 0, -1.29046011346], [str("4"), 0, 0, 0, 0, 1.71219173182, -1.61245472957], [ str("1_3"), 0, - 1.45632365305, + 1.45632365052, 0, - 0.821642099981, + 0.821642101161, 0, 0.0570637738537, ], @@ -826,9 +829,9 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( str("2_4"), 0, 0, - 1.79461739534, + 1.79461740051, 0, - 0.51578023092, + 0.515780222673, 0.433811455157, ], [ @@ -836,15 +839,15 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( 0, 0, 0, - 1.17191599967, - 0.93091304322, + 1.17191599323, + 0.930913047097, -0.739857285933, ], [ str("1_2"), 0, - 1.04119433365, - 1.44434211955, + 1.04119433442, + 1.444342119, 0, 0, 0.953575526751, @@ -852,19 +855,19 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( [ str("1_2_3"), 0, - 0.925047594605, - 1.32819621337, - 0.290365264803, + 0.925047594422, + 1.32819621359, + 0.290365264694, 0, 0.998573206243, ], [ str("1_2_3_4"), 0, - 0.914079355436, - 1.3172277155, - 0.279396891938, - 0.0383892587007, + 0.914079354784, + 1.31722771457, + 0.279396891914, + 0.0383892606964, 0.999295535579, ], ] @@ -920,8 +923,8 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_1"), - 2.0040856888, - -1.0020428444, + 2.00408568934, + -1.00204284467, 0, 0, 0, @@ -942,9 +945,9 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_2"), - 2.13847194795, + 2.13847194875, 0, - -1.06923597397, + -1.06923597438, 0, 0, 0, @@ -964,10 +967,10 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_3"), - 1.79252468631, + 1.79252468736, 0, 0, - -0.896262343155, + -0.89626234368, 0, 0, 0, @@ -1008,10 +1011,10 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_1_3"), - 2.27796575303, - -0.728161826525, + 2.27796575168, + -0.72816182526, 0, - -0.41082104999, + -0.41082105058, 0, 0, 0, @@ -1030,11 +1033,11 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_2_4"), - 2.31039762626, + 2.31039762318, 0, - -0.8973086976, + -0.897308700255, 0, - -0.25789011546, + -0.257890111337, 0, 0, 1, @@ -1052,11 +1055,11 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_3_4"), - 2.10282904289, + 2.10282904033, 0, 0, - -0.585957999835, - -0.46545652161, + -0.585957996615, + -0.46545652354, 0, 0, 1, @@ -1074,9 +1077,9 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_1_2"), - 2.4855364532, - -0.520597166825, - -0.722171059775, + 2.48553645342, + -0.52059716721, + -0.7221710595, 0, 0, 0, @@ -1096,10 +1099,10 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_1_2_3"), - 2.54360907278, - -0.462523797302, - -0.664098106685, - -0.145182632402, + 2.54360907271, + -0.462523797211, + -0.664098106795, + -0.145182632347, 0, 0, 0, @@ -1118,11 +1121,11 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH( ], [ str("0_1_2_3_4"), - 2.54909322157, - -0.457039677718, - -0.65861385775, - -0.139698445969, - -0.0191946293504, + 2.54909322196, + -0.457039677392, + -0.658613857285, + -0.139698445957, + -0.0191946303482, 0, 0, 1, @@ -1323,6 +1326,7 @@ def test_gaussian_style_files_fit_oplsaa_fit_CT_CT_C_OH_in_COOH_missing_1st_poin gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) assert ( @@ -1988,6 +1992,7 @@ def test_gaussian_style_files_fit_oplsaa_fit_CT_CT_C_OH_in_COOH(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=5e-03, + opls_force_k0_zero=True, ) assert ( @@ -2658,6 +2663,7 @@ def test_gaussian_style_files_fit_oplsaa_fit_CT_CT_C_OH_in_COOH_2_files_missing_ gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) assert ( @@ -3338,6 +3344,7 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH_bad gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_gaussian_style_files_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH_bad_element_order_mol2( @@ -3378,6 +3385,7 @@ def test_gaussian_style_files_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH_ gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=5e-03, + opls_force_k0_zero=True, ) def test_gaussian_log_file_variable_VDWGeometricSigma_default(self): @@ -3401,6 +3409,7 @@ def test_gaussian_log_file_variable_VDWGeometricSigma_default(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) with open( @@ -3443,6 +3452,7 @@ def test_gaussian_log_file_variable_VDWGeometricSigma_True(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) with open( @@ -3485,6 +3495,7 @@ def test_gaussian_log_file_variable_VDWGeometricSigma_False(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) with open( @@ -3533,6 +3544,35 @@ def test_bad_fit_dihedral_atom_types_input_list_of_3(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, + ) + + def test_bad_fit_opls_force_k0_zero_not_bool(self): + with pytest.raises( + TypeError, + match=r"ERROR: Please enter the 'opls_force_k0_zero' as a bool.", + ): + fit_dihedral_with_gomc( + ["HC", "CT", "CT"], + self.get_fn( + "gaussian/HC_CT_CT_HC/input/starting_coords/ethane_aa.mol2" + ), + self.get_fn("oplsaa_ethane_HC_CT_CT_HC.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian/HC_CT_CT_HC/output/HC_CT_CT_HC_multiplicity_1.log" + ): [0], + }, + zero_dihedral_atom_types=None, + qm_engine="gaussian", + combining_rule="lorentz", + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.99, + r_squared_rtol=1e-03, + opls_force_k0_zero="x", ) def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_0(self): @@ -3562,6 +3602,7 @@ def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_0(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_1(self): @@ -3591,6 +3632,7 @@ def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_1(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_2(self): @@ -3620,6 +3662,7 @@ def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_2(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_3(self): @@ -3649,6 +3692,7 @@ def test_bad_fit_dihedral_atom_types_input_list_of_4_with_int_at_3(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_mol2_file_file_does_not_exist(self): @@ -3676,6 +3720,7 @@ def test_mol2_file_file_does_not_exist(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_mol2_file_file_no_mol2_extention(self): @@ -3702,6 +3747,7 @@ def test_mol2_file_file_no_mol2_extention(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_mol2_file_file_not_a_string(self): @@ -3728,6 +3774,7 @@ def test_mol2_file_file_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_xml_selection_file_does_not_exist(self): @@ -3757,6 +3804,7 @@ def test_xml_selection_file_does_not_exist(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_xml_selection_file_no_xml_extention(self): @@ -3786,6 +3834,7 @@ def test_xml_selection_file_no_xml_extention(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_xml_selection_file_not_a_string(self): @@ -3814,6 +3863,7 @@ def test_xml_selection_file_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_temperature_unyt_units_not_a_temperture_but_pressure(self): @@ -3841,6 +3891,7 @@ def test_temperature_unyt_units_not_a_temperture_but_pressure(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_temperature_unyt_units_not_in_unyt_units(self): @@ -3868,6 +3919,7 @@ def test_temperature_unyt_units_not_in_unyt_units(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_not_a_dict(self): @@ -3893,6 +3945,7 @@ def test_qm_log_files_and_entries_not_a_dict(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_key_1_not_a_string(self): @@ -3923,6 +3976,7 @@ def test_qm_log_files_and_entries_key_1_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_key_2_not_a_string(self): @@ -3953,6 +4007,7 @@ def test_qm_log_files_and_entries_key_2_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_value_1_not_a_list(self): @@ -3985,6 +4040,7 @@ def test_qm_log_files_and_entries_value_1_not_a_list(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_value_2_not_a_list(self): @@ -4017,6 +4073,7 @@ def test_qm_log_files_and_entries_value_2_not_a_list(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_list_1_not_all_int(self): @@ -4049,6 +4106,7 @@ def test_qm_log_files_and_entries_list_1_not_all_int(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_list_2_not_all_int(self): @@ -4081,6 +4139,7 @@ def test_qm_log_files_and_entries_list_2_not_all_int(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_list_1_int_less_than_0(self): @@ -4113,6 +4172,7 @@ def test_qm_log_files_and_entries_list_1_int_less_than_0(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_qm_log_files_and_entries_list_2_int_less_than_0(self): @@ -4145,6 +4205,7 @@ def test_qm_log_files_and_entries_list_2_int_less_than_0(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_gomc_binary_path_not_a_string(self): @@ -4172,6 +4233,7 @@ def test_gomc_binary_path_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_gomc_binary_path_containing_the_GOMC_CPU_NVT_file_does_not_exist( @@ -4201,6 +4263,7 @@ def test_gomc_binary_path_containing_the_GOMC_CPU_NVT_file_does_not_exist( gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1e-03, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_not_list(self): @@ -4234,6 +4297,7 @@ def test_zero_dihedral_atom_types_not_list(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_1_str(self): @@ -4267,6 +4331,7 @@ def test_zero_dihedral_atom_types_list_1_str(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_2_str(self): @@ -4300,6 +4365,7 @@ def test_zero_dihedral_atom_types_list_2_str(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_1_not_4_strings(self): @@ -4336,6 +4402,7 @@ def test_zero_dihedral_atom_types_list_1_not_4_strings(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_2_not_4_strings(self): @@ -4372,6 +4439,7 @@ def test_zero_dihedral_atom_types_list_2_not_4_strings(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_1_not_lenght_4(self): @@ -4408,6 +4476,7 @@ def test_zero_dihedral_atom_types_list_1_not_lenght_4(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_zero_dihedral_atom_types_list_2_not_lenght_4(self): @@ -4444,6 +4513,7 @@ def test_zero_dihedral_atom_types_list_2_not_lenght_4(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_qm_engine_not_correct_value(self): @@ -4480,6 +4550,7 @@ def test_qm_engine_not_correct_value(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_qm_engine_not_a_string(self): @@ -4514,6 +4585,7 @@ def test_qm_engine_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_combining_rule_not_correct_value_is_list(self): @@ -4548,6 +4620,7 @@ def test_combining_rule_not_correct_value_is_list(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_combining_rule_not_correct_value_wrong_string(self): @@ -4582,6 +4655,7 @@ def test_combining_rule_not_correct_value_wrong_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_atom_type_naming_style_not_correct_value(self): @@ -4618,6 +4692,7 @@ def test_atom_type_naming_style_not_correct_value(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_atom_type_naming_style_not_a_string(self): @@ -4652,6 +4727,7 @@ def test_atom_type_naming_style_not_a_string(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_gomc_cpu_cores_not_correct_value(self): @@ -4686,6 +4762,7 @@ def test_gomc_cpu_cores_not_correct_value(self): gomc_cpu_cores=0, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_gomc_cpu_cores_not_a_int(self): @@ -4720,6 +4797,7 @@ def test_gomc_cpu_cores_not_a_int(self): gomc_cpu_cores=1.000, r_squared_min=0.99, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_r_squared_min_not_correct_value_is_0(self): @@ -4755,6 +4833,7 @@ def test_r_squared_min_not_correct_value_is_0(self): gomc_cpu_cores=1, r_squared_min=0.00, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_r_squared_min_not_correct_value_is_1(self): @@ -4790,6 +4869,7 @@ def test_r_squared_min_not_correct_value_is_1(self): gomc_cpu_cores=1, r_squared_min=1.00, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_r_squared_min_not_a_float(self): @@ -4825,6 +4905,7 @@ def test_r_squared_min_not_a_float(self): gomc_cpu_cores=1, r_squared_min=2, r_squared_rtol=0.02, + opls_force_k0_zero=True, ) def test_r_squared_rtol_not_correct_value_is_0(self): @@ -4860,6 +4941,7 @@ def test_r_squared_rtol_not_correct_value_is_0(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=0.00, + opls_force_k0_zero=True, ) def test_r_squared_rtol_not_correct_value_is_1(self): @@ -4895,6 +4977,7 @@ def test_r_squared_rtol_not_correct_value_is_1(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=1.00, + opls_force_k0_zero=True, ) def test_r_squared_rtol_not_a_float(self): @@ -4930,34 +5013,42 @@ def test_r_squared_rtol_not_a_float(self): gomc_cpu_cores=1, r_squared_min=0.99, r_squared_rtol=2, + opls_force_k0_zero=True, ) - def test_warning_r_squared_min_and_r_squared_rtol_need_adjusted( + def test_error_r_squared_min_and_r_squared_rtol_need_adjusted( self, ): with pytest.raises( ValueError, - match=f"ERROR: The calculated R-squared energy values from the fit type " - f"{'1_2_3'} " - f"does not match the validated case for 'r_squared_min' >= " + match=f"ERROR: The calculated R-squared energy values from the fit types " + f"do not match the validated case for 'r_squared_min' >= " f"{'0.99'}, " f"within the relative tolerance or 'r_squared_rtol' = " f"{'2e-07'}. \n" - f"- Fit via the individual or multi-dihedral fit, when " - f"Gaussian minus GOMC with the selected dihedral set to zero \n" - f"--> R-squared = " - f"{'0.998'} \n" - f"- Fit via the validation test case, when " - f"Gaussian minus GOMC with the selected individual dihedral added in GOMC \n" - f"-- >R-squared = " - f"{'0.987'} \n" + f"Constants_used = The calculated R-squared energy values from the fit type constants_used." + f"R-squared_fitted = Gaussian minus GOMC with the selected dihedral set to zero \n" + f"R-squared_new_dihedral = Gaussian minus GOMC with the selected individual dihedral added in GOMC \n" + f"Abs\(delta\) = absolute_value\(R-squared_fitted - R-squared_new_dihedral\) \n" + f"- The fits for all are shown here for all the dihedral combinations \n" + f"\[opls_constants_used, R-squared_fitted, R-squared_new_dihedral, Abs\(delta\)\] \n" + f"\[\['1', -0.27129793, -1.2824201, 1.0111222\], " + f"\['2', 0.45284933, -0.069821698, 0.52267103\], " + f"\['3', -1.5045949, -3.102543, 1.5979481\], " + f"\['4', -1.5448789, -3.1964352, 1.6515563\], " + f"\['1_3', -0.0060455503, -0.85802999, 0.85198444\], " + f"\['2_4', 0.55469005, 0.092830036, 0.46186001\], " + f"\['3_4', -0.8973394, -2.3260485, 1.4287091\], " + f"\['1_2', 0.98670161, 0.96338608, 0.02331553\], " + f"\['1_2_3', 0.99785712, 0.98674454, 0.01111258\], " + f"\['1_2_3_4', 0.99807697, 0.98725523, 0.01082174\]\], \n" f"The 'r_squared_min' and 'r_squared_rtol' " f"variables may need to be adjusted, \n" f"there is likely something wrong with the fitting procedure, the " f"software parameters need tuned, or there is a bug in the software. \n\n " f"NOTE: Since the R-squared values are calculated via different parameters, \n" f"the compared R-squared values could be very different if they are not nearly \n" - r"a perfect fit \(R-squared --> ~0.98 to 0.99999999\).", + f"a perfect fit \(R-squared --> ~0.98 to 0.99999999\).", ): fit_dihedral_with_gomc( ["CT", "CT", "C", "OH"], @@ -4985,6 +5076,57 @@ def test_warning_r_squared_min_and_r_squared_rtol_need_adjusted( r_squared_rtol=0.0000002, ) + def test_warning_r_squared_min_and_r_squared_rtol_need_adjusted( + self, + ): + with pytest.warns( + UserWarning, + match=f"WARNING: The calculated R-squared energy values from the fit types " + f"do match the validated case for 'r_squared_min' >= " + f"{'0.98'}, " + f"but do not fit within the relative tolerance of 'r_squared_rtol' = " + f"{'0.015'}. \n" + f"Constants_used = The calculated R-squared energy values from the fit type constants_used. \n" + f"R-squared_fitted = Gaussian minus GOMC with the selected dihedral set to zero. \n" + f"R-squared_new_dihedral = Gaussian minus GOMC with the selected individual dihedral added in GOMC. \n" + f"Abs\(delta\) = Abs\(R-squared_fitted - R-squared_new_dihedral\). \n" + f"- The fits for all are shown here for all the dihedral combinations \n" + f"\[opls_constants_used, R-squared_fitted, R-squared_new_dihedral, Abs\(delta\)\] \n" + f"\[\['1_2', 0.98670161, 0.96338608, 0.02331553\]\]. \n" + f"The 'r_squared_min' and 'r_squared_rtol' " + f"variables may need to be adjusted, \n" + f"there is likely something wrong with the fitting procedure, the " + f"software parameters need tuned, or there is a bug in the software. \n\n " + f"NOTE: Since the R-squared values are calculated via different parameters, \n" + f"the compared R-squared values could be very different if they are not nearly \n" + f"a perfect fit \(R-squared --> ~0.98 to 0.99999999\).", + ): + fit_dihedral_with_gomc( + ["CT", "CT", "C", "OH"], + self.get_fn( + "gaussian_style_output_files/CT_CT_C_OH/input/starting_coords/CT_CT_C_3_OH.mol2" + ), + self.get_fn("gmso_oplsaa_CT_CT_C_OH_in_COOH.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian_style_output_files/CT_CT_C_OH_split_part_2/output" + ): [], + self.get_fn( + "gaussian_style_output_files/CT_CT_C_OH_split_part_1/output" + ): [0], + }, + manual_dihedral_atom_numbers_list=[3, 2, 1, 4], + zero_dihedral_atom_types=[["CT", "CT", "C", "O_3"]], + qm_engine="gaussian_style_final_files", + combining_rule="geometric", + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.98, + r_squared_rtol=0.015, + ) + def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH_in_mie_form( self, ): @@ -5651,3 +5793,2686 @@ def test_gaussian_log_file_fit_oplsaa_protonated_fragment_CT_CT_C_OH_in_COOH_in_ atol=0.02, rtol=0.08, ) + + """ + """ + + def test_gaussian_style_files_fit_amber_aa_fit_CT_CT_CT_CT_in_butane_files( + self, + ): + fit_dihedral_with_gomc( + ["CT", "CT", "CT", "CT"], + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2" + ), + self.get_fn("amber_aa_butane_CT_CT_CT_CT_gmso.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/output" + ): [], + }, + manual_dihedral_atom_numbers_list=[1, 2, 3, 4], + zero_dihedral_atom_types=None, + qm_engine="gaussian_style_final_files", + combining_rule=None, + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.98, + r_squared_rtol=0.02, + opls_force_k0_zero=False, + ) + + assert ( + os.path.isfile("all_normalized_energies_in_kcal_per_mol.txt") + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_4_in_kcal_per_mol.txt" + ) + is True + ) + assert os.path.isfile("gomc_raw_energies_in_Kelvin.txt") is True + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_4_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile( + "opls_all_single_fit_dihedral_k_constants_figure.pdf" + ) + is True + ) + assert ( + os.path.isfile("opls_all_summed_dihedrals_k_constants_figure.pdf") + is True + ) + assert ( + os.path.isfile("opls_dihedral_k_constants_fit_energy.txt") is True + ) + assert ( + os.path.isfile("periodic_dihedral_k_constants_fit_energy.txt") + is True + ) + assert os.path.isfile("RB_torsion_k_constants_fit_energy.txt") is True + + # check the OPLS dihedral file + with open("opls_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "r_squared", + ], + [str("1"), 0, 2.13492829593, 0, 0, 0, -0.851303269629], + [str("2"), 0, 0, 2.25256162343, 0, 0, -0.690454186065], + [ + str("3"), + 6.12246736014, + 0, + 0, + -3.06123368007, + 0, + 0.960991418867, + ], + [str("4"), 0, 0, 0, 0, 2.0884153941, -0.91251897505], + [ + str("1_3"), + 6.5396241564, + -0.208579929672, + 0, + -3.06123214853, + 0, + 0.977746134654, + ], + [ + str("2_4"), + 0, + 0, + 1.54851017553, + 0, + 1.05607698155, + -0.497352492204, + ], + [ + str("3_4"), + 6.66816846056, + 0, + 0, + -3.06123287378, + -0.348137574352, + 0.985623284101, + ], + [ + str("1_2"), + 0, + 1.13979595936, + 1.49269914954, + 0, + 0, + -0.465523847167, + ], + [ + str("1_2_3"), + 6.53963027886, + -0.208580993867, + 0.144337057476, + -3.06123414556, + 0, + 0.977793543598, + ], + [ + str("1_2_3_4"), + 6.55860366135, + -0.208578022741, + 0.14433854159, + -3.06123334065, + -0.348136431047, + 0.982572116976, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the periodic dihedral file + with open("periodic_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "n0_kcal_per_mol", + "n1_kcal_per_mol", + "n2_kcal_per_mol", + "n3_kcal_per_mol", + "n4_kcal_per_mol", + "n5_kcal_per_mol", + "d0_kcal_per_mol", + "d1_kcal_per_mol", + "d2_kcal_per_mol", + "d3_kcal_per_mol", + "d4_kcal_per_mol", + "d5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 2.1349282959, + -1.06746414796, + 0, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.851303269629, + ], + [ + str("0_2"), + 2.25256162343, + 0, + -1.12628081172, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.690454186065, + ], + [ + str("0_3"), + 0, + 0, + 0, + 1.53061684004, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.960991418867, + ], + [ + str("0_4"), + 2.0884153941, + 0, + 0, + 0, + -1.04420769705, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.91251897505, + ], + [ + str("0_1_3"), + -1.99973371195e-12, + 0.104289964836, + 0, + 1.53061607426, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.977746134654, + ], + [ + str("0_2_4"), + 2.60458715708, + 0, + -0.774255087765, + 0, + -0.528038490775, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.497352492204, + ], + [ + str("0_3_4"), + -0.075286217852, + 0, + 0, + 1.53061643689, + 0.174068787176, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.985623284101, + ], + [ + str("0_1_2"), + 2.6324951089, + -0.56989797968, + -0.74634957477, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.465523847167, + ], + [ + str("0_1_2_3"), + 0.144337057479, + 0.104290496934, + -0.072168528738, + 1.53061707278, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.977793543598, + ], + [ + str("0_1_2_3_4"), + -0.194307422173, + 0.10428901137, + -0.072169270795, + 1.53061667032, + 0.174068215524, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.982572116976, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the RB torsion file + with open("RB_torsion_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 1.06746414796, + -1.06746414796, + 0, + 0, + 0, + 0, + -0.851303269629, + ], + [ + str("0_2"), + 2.25256162343, + 0, + -2.25256162343, + 0, + 0, + 0, + -0.690454186065, + ], + [ + str("0_3"), + 1.53061684004, + -4.5918505201, + 0, + 6.12246736014, + 0, + 0, + 0.960991418867, + ], + [ + str("0_4"), + 0, + 0, + 8.3536615764, + 0, + -8.3536615764, + 0, + -0.91251897505, + ], + [ + str("0_1_3"), + 1.6349060391, + -4.48755825796, + 0, + 6.12246429706, + 0, + 0, + 0.977746134654, + ], + [ + str("0_2_4"), + 1.54851017553, + 0, + 2.67579775067, + 0, + -4.2243079262, + 0, + -0.497352492204, + ], + [ + str("0_3_4"), + 1.80346779339, + -4.59184931067, + -1.39255029741, + 6.12246574756, + 1.39255029741, + 0, + 0.985623284101, + ], + [ + str("0_1_2"), + 2.06259712922, + -0.56989797968, + -1.49269914954, + 0, + 0, + 0, + -0.465523847167, + ], + [ + str("0_1_2_3"), + 1.77924462719, + -4.48756072141, + -0.144337057476, + 6.12246829112, + 0, + 0, + 0.977793543598, + ], + [ + str("0_1_2_3_4"), + 1.78873469057, + -4.4875609996, + -1.53688426578, + 6.1224666813, + 1.39254572419, + 0, + 0.982572116976, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + def test_gaussian_style_files_fit_opls_ua_fit_CT_CT_CT_CT_in_butane_files( + self, + ): + fit_dihedral_with_gomc( + ["CH3", "CH2", "CH2", "CH3"], + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2" + ), + self.get_fn("trappe_ua_butane_CT_CT_CT_CT_gmso.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/output" + ): [], + }, + manual_dihedral_atom_numbers_list=[1, 2, 3, 4], + zero_dihedral_atom_types=None, + qm_engine="gaussian_style_final_files", + combining_rule=None, + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.97, + r_squared_rtol=0.02, + opls_force_k0_zero=True, + ) + + assert ( + os.path.isfile("all_normalized_energies_in_kcal_per_mol.txt") + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_4_in_kcal_per_mol.txt" + ) + is True + ) + assert os.path.isfile("gomc_raw_energies_in_Kelvin.txt") is True + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_4_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile( + "opls_all_single_fit_dihedral_k_constants_figure.pdf" + ) + is True + ) + assert ( + os.path.isfile("opls_all_summed_dihedrals_k_constants_figure.pdf") + is True + ) + assert ( + os.path.isfile("opls_dihedral_k_constants_fit_energy.txt") is True + ) + assert ( + os.path.isfile("periodic_dihedral_k_constants_fit_energy.txt") + is True + ) + assert os.path.isfile("RB_torsion_k_constants_fit_energy.txt") is True + + # check the OPLS dihedral file + with open("opls_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "r_squared", + ], + [str("1"), 0, 3.89019135095, 0, 0, 0, -0.15385629089], + [str("2"), 0, 0, 3.03580389236, 0, 0, -1.08874124105], + [str("3"), 0, 0, 0, 4.57026969264, 0, 0.755185044925], + [str("4"), 0, 0, 0, 0, 3.24432472808, -0.881847591138], + [ + str("1_3"), + 0, + 1.51800454722, + 0, + 3.55826890289, + 0, + 0.957438274954, + ], + [ + str("2_4"), + 0, + 0, + 1.57125480195, + 0, + 2.19682320759, + -0.665155000765, + ], + [ + str("3_4"), + 0, + 0, + 0, + 4.33329974254, + 0.355455279474, + 0.76627477714, + ], + [ + str("1_2"), + 0, + 3.35937487818, + 0.796225057519, + 0, + 0, + -0.098211712848, + ], + [ + str("1_2_3"), + 0, + 1.81662366124, + -0.74655694329, + 3.85689590923, + 0, + 0.998529953567, + ], + [ + str("1_2_3_4"), + 0, + 1.84350327627, + -0.719677798851, + 3.88377504842, + -0.0940775749305, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the periodic dihedral file + with open("periodic_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "n0_kcal_per_mol", + "n1_kcal_per_mol", + "n2_kcal_per_mol", + "n3_kcal_per_mol", + "n4_kcal_per_mol", + "n5_kcal_per_mol", + "d0_kcal_per_mol", + "d1_kcal_per_mol", + "d2_kcal_per_mol", + "d3_kcal_per_mol", + "d4_kcal_per_mol", + "d5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 3.89019135095, + -1.94509567548, + 0, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.15385629089, + ], + [ + str("0_2"), + 3.03580389236, + 0, + -1.51790194618, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -1.08874124105, + ], + [ + str("0_3"), + 4.57026969264, + 0, + 0, + -2.28513484632, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.755185044925, + ], + [ + str("0_4"), + 3.24432472808, + 0, + 0, + 0, + -1.62216236404, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.881847591138, + ], + [ + str("0_1_3"), + 5.07627345011, + -0.75900227361, + 0, + -1.77913445144, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.957438274954, + ], + [ + str("0_2_4"), + 3.76807800954, + 0, + -0.785627400975, + 0, + -1.0984116038, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.665155000765, + ], + [ + str("0_3_4"), + 4.68875502201, + 0, + 0, + -2.16664987127, + -0.177727639737, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.76627477714, + ], + [ + str("0_1_2"), + 4.1555999357, + -1.67968743909, + -0.39811252876, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.098211712848, + ], + [ + str("0_1_2_3"), + 4.92696262718, + -0.90831183062, + 0.373278471645, + -1.92844795462, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.998529953567, + ], + [ + str("0_1_2_3_4"), + 4.91352295091, + -0.921751638135, + 0.359838899426, + -1.94188752421, + 0.0470387874652, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the RB torsion file + with open("RB_torsion_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 1.94509567548, + -1.94509567548, + 0, + 0, + 0, + 0, + -0.15385629089, + ], + [ + str("0_2"), + 3.03580389236, + 0, + -3.03580389236, + 0, + 0, + 0, + -1.08874124105, + ], + [ + str("0_3"), + 2.28513484632, + 6.85540453896, + 0, + -9.14053938528, + 0, + 0, + 0.755185044925, + ], + [ + str("0_4"), + 0, + 0, + 12.9772989123, + 0, + -12.9772989123, + 0, + -0.881847591138, + ], + [ + str("0_1_3"), + 2.53813672505, + 4.57840108073, + 0, + -7.11653780578, + 0, + 0, + 0.957438274954, + ], + [ + str("0_2_4"), + 1.57125480195, + 0, + 7.21603802841, + 0, + -8.78729283036, + 0, + -0.665155000765, + ], + [ + str("0_3_4"), + 2.16664987127, + 6.49994961381, + 1.4218211179, + -8.66659948508, + -1.4218211179, + 0, + 0.76627477714, + ], + [ + str("0_1_2"), + 2.47591249661, + -1.67968743909, + -0.796225057519, + 0, + 0, + 0, + -0.098211712848, + ], + [ + str("0_1_2_3"), + 2.09020284194, + 4.877032033237, + 0.74655694329, + -7.71379181846, + 0, + 0, + 0.998529953567, + ], + [ + str("0_1_2_3_4"), + 2.14396136349, + 4.90391093449, + 0.343367499129, + -7.76755009684, + 0.376310299722, + 0, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + def test_gaussian_style_files_fit_mia_ua_fit_CT_CT_CT_CT_in_butane_files( + self, + ): + fit_dihedral_with_gomc( + ["CH3", "CH2", "CH2", "CH3"], + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2" + ), + self.get_fn("mie_ua_butane_CT_CT_CT_CT_gmso.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/output" + ): [], + }, + manual_dihedral_atom_numbers_list=[1, 2, 3, 4], + zero_dihedral_atom_types=None, + qm_engine="gaussian_style_final_files", + combining_rule=None, + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.98, + r_squared_rtol=0.02, + opls_force_k0_zero=False, + ) + + assert ( + os.path.isfile("all_normalized_energies_in_kcal_per_mol.txt") + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_4_in_kcal_per_mol.txt" + ) + is True + ) + assert os.path.isfile("gomc_raw_energies_in_Kelvin.txt") is True + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_4_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile( + "opls_all_single_fit_dihedral_k_constants_figure.pdf" + ) + is True + ) + assert ( + os.path.isfile("opls_all_summed_dihedrals_k_constants_figure.pdf") + is True + ) + assert ( + os.path.isfile("opls_dihedral_k_constants_fit_energy.txt") is True + ) + assert ( + os.path.isfile("periodic_dihedral_k_constants_fit_energy.txt") + is True + ) + assert os.path.isfile("RB_torsion_k_constants_fit_energy.txt") is True + + # check the OPLS dihedral file + with open("opls_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "r_squared", + ], + [str("1"), 0, 3.89019135095, 0, 0, 0, -0.15385629089], + [str("2"), 0, 0, 3.03580389236, 0, 0, -1.08874124105], + [str("3"), 0, 0, 0, 4.57026969264, 0, 0.755185044925], + [str("4"), 0, 0, 0, 0, 3.24432472808, -0.881847591138], + [ + str("1_3"), + 0, + 1.51800454722, + 0, + 3.55826890289, + 0, + 0.957438274954, + ], + [ + str("2_4"), + 0, + 0, + 1.57125480195, + 0, + 2.19682320759, + -0.665155000765, + ], + [ + str("3_4"), + 0, + 0, + 0, + 4.33329974254, + 0.355455279474, + 0.76627477714, + ], + [ + str("1_2"), + 0, + 3.35937487818, + 0.796225057519, + 0, + 0, + -0.098211712848, + ], + [ + str("1_2_3"), + 0, + 1.81662366124, + -0.74655694329, + 3.85689590923, + 0, + 0.998529953567, + ], + [ + str("1_2_3_4"), + 0, + 1.84350327627, + -0.719677798851, + 3.88377504842, + -0.0940775749305, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the periodic dihedral file + with open("periodic_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "n0_kcal_per_mol", + "n1_kcal_per_mol", + "n2_kcal_per_mol", + "n3_kcal_per_mol", + "n4_kcal_per_mol", + "n5_kcal_per_mol", + "d0_kcal_per_mol", + "d1_kcal_per_mol", + "d2_kcal_per_mol", + "d3_kcal_per_mol", + "d4_kcal_per_mol", + "d5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 3.89019135095, + -1.94509567548, + 0, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.15385629089, + ], + [ + str("0_2"), + 3.03580389236, + 0, + -1.51790194618, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -1.08874124105, + ], + [ + str("0_3"), + 4.57026969264, + 0, + 0, + -2.28513484632, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.755185044925, + ], + [ + str("0_4"), + 3.24432472808, + 0, + 0, + 0, + -1.62216236404, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.881847591138, + ], + [ + str("0_1_3"), + 5.07627345011, + -0.75900227361, + 0, + -1.77913445144, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.957438274954, + ], + [ + str("0_2_4"), + 3.76807800954, + 0, + -0.785627400975, + 0, + -1.0984116038, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.665155000765, + ], + [ + str("0_3_4"), + 4.68875502201, + 0, + 0, + -2.16664987127, + -0.177727639737, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.76627477714, + ], + [ + str("0_1_2"), + 4.1555999357, + -1.67968743909, + -0.39811252876, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.098211712848, + ], + [ + str("0_1_2_3"), + 4.92696262718, + -0.90831183062, + 0.373278471645, + -1.92844795462, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.998529953567, + ], + [ + str("0_1_2_3_4"), + 4.91352295091, + -0.921751638135, + 0.359838899426, + -1.94188752421, + 0.0470387874652, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the RB torsion file + with open("RB_torsion_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 1.94509567548, + -1.94509567548, + 0, + 0, + 0, + 0, + -0.15385629089, + ], + [ + str("0_2"), + 3.03580389236, + 0, + -3.03580389236, + 0, + 0, + 0, + -1.08874124105, + ], + [ + str("0_3"), + 2.28513484632, + 6.85540453896, + 0, + -9.14053938528, + 0, + 0, + 0.755185044925, + ], + [ + str("0_4"), + 0, + 0, + 12.9772989123, + 0, + -12.9772989123, + 0, + -0.881847591138, + ], + [ + str("0_1_3"), + 2.53813672505, + 4.57840108073, + 0, + -7.11653780578, + 0, + 0, + 0.957438274954, + ], + [ + str("0_2_4"), + 1.57125480195, + 0, + 7.21603802841, + 0, + -8.78729283036, + 0, + -0.665155000765, + ], + [ + str("0_3_4"), + 2.16664987127, + 6.49994961381, + 1.4218211179, + -8.66659948508, + -1.4218211179, + 0, + 0.76627477714, + ], + [ + str("0_1_2"), + 2.47591249661, + -1.67968743909, + -0.796225057519, + 0, + 0, + 0, + -0.098211712848, + ], + [ + str("0_1_2_3"), + 2.09020284194, + 4.877032033237, + 0.74655694329, + -7.71379181846, + 0, + 0, + 0.998529953567, + ], + [ + str("0_1_2_3_4"), + 2.14396136349, + 4.90391093449, + 0.343367499129, + -7.76755009684, + 0.376310299722, + 0, + 0.999129219573, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + def test_gaussian_style_files_fit_amber_fit_CT_CT_CT_CT_in_butane_files( + self, + ): + fit_dihedral_with_gomc( + ["CT", "CT", "CT", "CT"], + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/input/starting_coords/butane_aa.mol2" + ), + self.get_fn("exp6_aa_butane_CT_CT_CT_CT_gmso.xml"), + 298.15 * u.Kelvin, + gomc_binary_directory, + { + self.get_fn( + "gaussian_style_output_files/CT_CT_CT_CT/output" + ): [], + }, + manual_dihedral_atom_numbers_list=[1, 2, 3, 4], + zero_dihedral_atom_types=None, + qm_engine="gaussian_style_final_files", + combining_rule=None, + atom_type_naming_style="general", + gomc_cpu_cores=1, + r_squared_min=0.98, + r_squared_rtol=0.02, + opls_force_k0_zero=False, + ) + + assert ( + os.path.isfile("all_normalized_energies_in_kcal_per_mol.txt") + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_1_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_2_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_4_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_3_in_kcal_per_mol.txt" + ) + is True + ) + assert ( + os.path.isfile( + "all_normalized_energies_OPLS_fit_4_in_kcal_per_mol.txt" + ) + is True + ) + assert os.path.isfile("gomc_raw_energies_in_Kelvin.txt") is True + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_2_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_3_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_1_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_2_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_4_energies_in_Kelvin.txt") + is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_3_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile("gomc_raw_OPLS_fit_4_energies_in_Kelvin.txt") is True + ) + assert ( + os.path.isfile( + "opls_all_single_fit_dihedral_k_constants_figure.pdf" + ) + is True + ) + assert ( + os.path.isfile("opls_all_summed_dihedrals_k_constants_figure.pdf") + is True + ) + assert ( + os.path.isfile("opls_dihedral_k_constants_fit_energy.txt") is True + ) + assert ( + os.path.isfile("periodic_dihedral_k_constants_fit_energy.txt") + is True + ) + assert os.path.isfile("RB_torsion_k_constants_fit_energy.txt") is True + + # check the OPLS dihedral file + with open("opls_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "r_squared", + ], + [str("1"), 0, 4.64189989445, 0, 0, 0, 0.0554910823082], + [str("2"), 0, 0, 3.62275916028, 0, 0, -1.23591651342], + [str("3"), 0, 0, 0, 4.92948888627, 0, 0.477550076455], + [str("4"), 0, 0, 0, 0, 3.7172458491, -1.12957843414], + [ + str("1_3"), + 0, + 2.4400162262, + 0, + 3.30281500888, + 0, + 0.9846783871, + ], + [ + str("2_4"), + 0, + 0, + 2.06026818365, + 0, + 2.34373594682, + -0.768017540175, + ], + [ + str("3_4"), + 0, + 0, + 0, + 4.41238953459, + 0.775649812101, + 0.52879674349, + ], + [ + str("1_2"), + 0, + 4.00810304165, + 0.950695705123, + 0, + 0, + 0.132478073226, + ], + [ + str("1_2_3"), + 0, + 2.61642006674, + -0.441015014606, + 3.47922351339, + 0, + 0.998594497652, + ], + [ + str("1_2_3_4"), + 0, + 2.65142528583, + -0.406010416441, + 3.51422811225, + -0.122516855939, + 0.999580827406, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the periodic dihedral file + with open("periodic_dihedral_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "n0_kcal_per_mol", + "n1_kcal_per_mol", + "n2_kcal_per_mol", + "n3_kcal_per_mol", + "n4_kcal_per_mol", + "n5_kcal_per_mol", + "d0_kcal_per_mol", + "d1_kcal_per_mol", + "d2_kcal_per_mol", + "d3_kcal_per_mol", + "d4_kcal_per_mol", + "d5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 4.64189989445, + -2.32094994722, + 0, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.0554910823082, + ], + [ + str("0_2"), + 3.62275916028, + 0, + -1.81137958014, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -1.23591651342, + ], + [ + str("0_3"), + 4.92948888627, + 0, + 0, + -2.46474444314, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.477550076455, + ], + [ + str("0_4"), + 3.7172458491, + 0, + 0, + 0, + -1.85862292455, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -1.12957843414, + ], + [ + str("0_1_3"), + 5.74283123508, + -1.2200081131, + 0, + -1.65140750444, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.9846783871, + ], + [ + str("0_2_4"), + 4.40400413047, + 0, + -1.03013409182, + 0, + -1.17186797341, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + -0.768017540175, + ], + [ + str("0_3_4"), + 5.18803934669, + 0, + 0, + -2.2061947673, + -0.38782490605, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.52879674349, + ], + [ + str("0_1_2"), + 4.95879874677, + -2.00405152082, + -0.475347852562, + 0, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.132478073226, + ], + [ + str("0_1_2_3"), + 5.65462856552, + -1.3082100333, + 0.220507507303, + -1.7396117567, + 0, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.998594497652, + ], + [ + str("0_1_2_3_4"), + 5.6371261257, + -1.32571264291, + 0.20300520822, + -1.75711405612, + 0.0612584279695, + 0, + 0, + 1, + 2, + 3, + 4, + 5, + 90.0, + 180.0, + 0, + 180.0, + 0, + 180.0, + 0.999580827406, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) + + # check the RB torsion file + with open("RB_torsion_k_constants_fit_energy.txt", "r") as fp: + number_sig_i = 4 + out_gomc = fp.readlines() + for i, line in enumerate(out_gomc): + split_line_i = line.split() + correct_line_values = [ + [ + "non_zero_k_constants", + "k0_kcal_per_mol", + "k1_kcal_per_mol", + "k2_kcal_per_mol", + "k3_kcal_per_mol", + "k4_kcal_per_mol", + "k5_kcal_per_mol", + "r_squared", + ], + [ + str("0_1"), + 2.32094994722, + -2.32094994722, + 0, + 0, + 0, + 0, + 0.0554910823082, + ], + [ + str("0_2"), + 3.62275916028, + 0, + -3.62275916028, + 0, + 0, + 0, + -1.23591651342, + ], + [ + str("0_3"), + 2.46474444314, + 7.3942333294, + 0, + -9.85897777254, + 0, + 0, + 0.477550076455, + ], + [ + str("0_4"), + 0, + 0, + 14.8689833964, + 0, + -14.8689833964, + 0, + -1.12957843414, + ], + [ + str("0_1_3"), + 2.87141561754, + 3.73421440022, + 0, + -6.60563001776, + 0, + 0, + 0.9846783871, + ], + [ + str("0_2_4"), + 2.06026818365, + 0, + 7.31467560363, + 0, + -9.37494378728, + 0, + -0.768017540175, + ], + [ + str("0_3_4"), + 2.2061947673, + 6.61858430188, + 3.1025992484, + -8.82477906918, + -3.1025992484, + 0, + 0.52879674349, + ], + [ + str("0_1_2"), + 2.95474722595, + -2.00405152082, + -0.950695705123, + 0, + 0, + 0, + 0.132478073226, + ], + [ + str("0_1_2_3"), + 2.60680677546, + 3.91062523672, + 0.441015014606, + -6.95844702678, + 0, + 0, + 0.998594497652, + ], + [ + str("0_1_2_3_4"), + 2.6768162826, + 3.94562952546, + -0.084057007315, + -7.0284562245, + 0.490067423756, + 0, + 0.999580827406, + ], + ] + + if i == 0: + assert split_line_i == correct_line_values[i] + + else: + assert len(split_line_i) == len(correct_line_values[i]) + for j in range(0, len(correct_line_values[i])): + # check the string listing what cosin powers are used + if j == 0: + assert str(split_line_i[j]) == str( + correct_line_values[i][j] + ) + + # check the k-values and the r-squared fit + else: + assert np.isclose( + mdf_math.round_to_sig_figs( + float(split_line_i[j]), number_sig_i + ), + mdf_math.round_to_sig_figs( + correct_line_values[i][j], number_sig_i + ), + atol=0.02, + rtol=0.08, + ) diff --git a/mosdef_dihedral_fit/utils/math_operations.py b/mosdef_dihedral_fit/utils/math_operations.py index 678ce3d..8956f81 100755 --- a/mosdef_dihedral_fit/utils/math_operations.py +++ b/mosdef_dihedral_fit/utils/math_operations.py @@ -275,7 +275,7 @@ def sum_opls_const_1_plus_or_minus_cos_n_values(phi_list): that the dihedral is fit too. Example 0: - const_1_minus_Cos_0_phi = 0 , since k0 is not used in this form + const_1_minus_Cos_0_phi = 1 , since k0 is a constant Example 1: const_1_plus_Cos_1_phi = sum of the list using all phis [(1 + cos(1 * phi))] in k1 * (1 + cos(1 * phi)) @@ -296,6 +296,7 @@ def sum_opls_const_1_plus_or_minus_cos_n_values(phi_list): Standard OPLS dihedral form .. math:: opls_dihedral &= 1/2 *( + &= k0 &= + k1 * (1 + cos(1 * phi)) &= + k2 * (1 - cos(2 * phi)) &= + k3 * (1 + cos(3 * phi)) @@ -304,6 +305,7 @@ def sum_opls_const_1_plus_or_minus_cos_n_values(phi_list): Modified OPLS dihedral form for all dihedrals with the same atom types/classes exist in the molecule. opls_dihedral_n_1 &= 1/2 *( + &= k0 &= + k1 * const_1_plus_Cos_1_phi &= + k2 * const_1_minus_Cos_2_phi &= + k3 * const_1_plus_Cos_3_phi @@ -329,7 +331,7 @@ def sum_opls_const_1_plus_or_minus_cos_n_values(phi_list): const_1_minus_Cos_4_phi: float (unitless) The sum of all phi values in the list [(1 - cos(4 * phi))] in k4 * (1 - cos(4 * phi)) """ - const_1_minus_Cos_0_phi = 0 + const_1_minus_Cos_0_phi = 1 const_1_plus_Cos_1_phi = 0 const_1_minus_Cos_2_phi = 0 const_1_plus_Cos_3_phi = 0 @@ -366,8 +368,6 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): The k-values are in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...), and the output ' dihedral_energy' energy are output in the same energy units. - NOTE: In this case, k0 is always set to zero (0) regardless of the - user entered value, because it is not in the equation form. NOTE: THIS WILL FIT THE FORM IT IS FED, REGARDLESS IF IT IS THE CORRECT ANALYTICAL SOLUTION. MEANING SOMETIMES YOU CAN NOT USE A K-VALUE BECAUSE @@ -378,6 +378,7 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): .. math:: opls_dihedral_n_1 &= 1/2 *( + &= k0 &= + k1 * (1 + cos(1 * phi)) &= + k2 * (1 - cos(2 * phi)) &= + k3 * (1 + cos(3 * phi)) @@ -386,7 +387,7 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): Parameters ---------- - cos_powers_phi_and_constants_data: tuple of (cos_powers, scanned_phi, all_sum_opls_const_1_plus_or_minus_cos_n_list) + cos_powers_phi_and_constants_data: list or tuple of (cos_powers, scanned_phi, all_sum_opls_const_1_plus_or_minus_cos_n_list) cos_powers: str, int, or float, The str options are: '1', '1_3', '2', '2_4', '1_2', '1_2_3', and '1_2_3_4' The int options are: '1', '13', '2', '24', '12', '123', and '1234' @@ -402,7 +403,7 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): Example 5: '1' only set the k1 to a non-zero value Example 6: '123.' only set the k1, k2, and k3 to a non-zero value - scanned_phi: floats, int, nest list of floats/int + scanned_phi: floats, int, nest list or tuple of floats/int The 'phi_data' angle of the dihedral is in degrees The floats, int options: the phi_data from a single point @@ -411,15 +412,7 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): being rotated in QM, and all the other identical dihedral angles from the same atom typed/classed atoms in the test molecule/fragment. - const_1_minus_Cos_0_phi: list, default=None - A list of the OPLS k0 data, which is always zero in this case - - Example 0: - const_1_minus_Cos_0_phi = 0 , since k0 is not used in this form - - const_1_plus_Cos_1_phi: values for all k-values, which - is required to fit the data, especially when multiple dihedrals of the same - atom types/classes exist in the molecule that the dihedral is fit too. + all_sum_opls_const_1_plus_or_minus_cos_n_list: list or tuple, default=None all_sum_opls_const_1_plus_or_minus_cos_n_list = [ const_1_minus_Cos_0_phi, @@ -429,8 +422,19 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): const_1_minus_Cos_4_phi ] + A list of the OPLS k0 data, which is always 0 or 1 in this case. + If all 0 values --> (0, 0, .., 0), then k0=0. + If all 1 values --> (1, 1, .., 1), then k0=constant + + # It is critical that 'const_1_minus_Cos_0_phi' be all (0, 0, .., 0) for k0=0 + # and all (1, 1, .., 1) 'const_1_minus_Cos_0_phi' for k0=constant + + const_1_plus_Cos_1_phi: values for all k-values, which + is required to fit the data, especially when multiple dihedrals of the same + atom types/classes exist in the molecule that the dihedral is fit too. + Example 0: - const_1_minus_Cos_0_phi = 0 , since k0 is not used in this form + const_1_minus_Cos_0_phi = k0 Example 1: const_1_plus_Cos_1_phi = sum of all phi values in the list @@ -454,9 +458,10 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): k0: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) The 'k0' value is the k-value for the opls dihedral where n=0, - or the constant without a cosine multiple. - NOTE: In this case, it is set to zero (0) regardless of the - user entered value, because it is not in the equation form. + or the constant without a cosine multiple. The 'k0' can be a + constant or zero (0). + Please see 'all_sum_opls_const_1_plus_or_minus_cos_n_list' variable + for more details. k1: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) The 'k1' value is the k-value for the opls dihedral where n=1 in the cosine multiple. @@ -536,14 +541,12 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): if ( const_1_minus_Cos_0_phi_data is None - or const_1_minus_Cos_0_phi_data is None or const_1_plus_Cos_1_phi_data is None or const_1_minus_Cos_2_phi_data is None or const_1_plus_Cos_3_phi_data is None or const_1_minus_Cos_4_phi_data is None ) and ( const_1_minus_Cos_0_phi_data - != const_1_minus_Cos_0_phi_data != const_1_plus_Cos_1_phi_data != const_1_minus_Cos_2_phi_data != const_1_plus_Cos_3_phi_data @@ -676,40 +679,55 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): f" 1., 2., 3., 4., 13., 24., 12., 34., 123., 1234.]" ) + # check if using k0=0 ('opls_force_k0_zero'=True) by summing + # 'sorted_const_1_minus_Cos_0_phi_data_lists' + # It is critical that 'const_1_minus_Cos_0_phi_data' be all (0, 0, .., 0) for k0=0 + # and all (1, 1, .., 1) 'const_1_minus_Cos_0_phi_data' for k0=constant + if ( + const_1_minus_Cos_0_phi_data is not None + and sum(const_1_minus_Cos_0_phi_data) == 0 + # or const_1_minus_Cos_0_phi_data is None + ): + k0 = 0 + if cos_powers_modified in ["1", 1, 1.0]: # NOTE: THE ALL BUT THE 'kx' VALUES ARE REPLACES WITH A CONSTANT, if use_const_1_plus_minus_Cos_x_values_bool is False: dihedral_energy = ( - 1 / 2 * (+k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180))) + 1 / 2 * (k0 + k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180))) ) elif use_const_1_plus_minus_Cos_x_values_bool is True: - dihedral_energy = 1 / 2 * (+k1 * const_1_plus_Cos_1_phi_data) + dihedral_energy = 1 / 2 * (k0 + k1 * const_1_plus_Cos_1_phi_data) elif cos_powers_modified in ["2", 2, 2.0]: if use_const_1_plus_minus_Cos_x_values_bool is False: dihedral_energy = ( - 1 / 2 * (+k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180))) + 1 / 2 * (k0 + k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180))) ) + elif use_const_1_plus_minus_Cos_x_values_bool is True: - dihedral_energy = 1 / 2 * (+k2 * const_1_minus_Cos_2_phi_data) + dihedral_energy = 1 / 2 * (k0 + k2 * const_1_minus_Cos_2_phi_data) elif cos_powers_modified in ["3", 3, 3.0]: if use_const_1_plus_minus_Cos_x_values_bool is False: dihedral_energy = ( - 1 / 2 * (+k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180))) + 1 / 2 * (k0 + k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180))) ) + elif use_const_1_plus_minus_Cos_x_values_bool is True: - dihedral_energy = 1 / 2 * (+k3 * const_1_plus_Cos_3_phi_data) + + dihedral_energy = 1 / 2 * (k0 + k3 * const_1_plus_Cos_3_phi_data) elif cos_powers_modified in ["4", 4, 4.0]: if use_const_1_plus_minus_Cos_x_values_bool is False: dihedral_energy = ( - 1 / 2 * (+k4 * (1 - np.cos(4 * scanned_phi * np.pi / 180))) + 1 / 2 * (k0 + k4 * (1 - np.cos(4 * scanned_phi * np.pi / 180))) ) + elif use_const_1_plus_minus_Cos_x_values_bool is True: - dihedral_energy = 1 / 2 * (+k4 * const_1_minus_Cos_4_phi_data) + dihedral_energy = 1 / 2 * (k0 + k4 * const_1_minus_Cos_4_phi_data) elif cos_powers_modified in ["1_3", 13, 13.0]: if use_const_1_plus_minus_Cos_x_values_bool is False: @@ -717,16 +735,19 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k0 + + k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180)) ) ) + elif use_const_1_plus_minus_Cos_x_values_bool is True: dihedral_energy = ( 1 / 2 * ( - +k1 * const_1_plus_Cos_1_phi_data + k0 + + +k1 * const_1_plus_Cos_1_phi_data + k3 * const_1_plus_Cos_3_phi_data ) ) @@ -737,7 +758,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180)) + k0 + + +k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180)) + k4 * (1 - np.cos(4 * scanned_phi * np.pi / 180)) ) ) @@ -747,7 +769,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k2 * const_1_minus_Cos_2_phi_data + k0 + + +k2 * const_1_minus_Cos_2_phi_data + k4 * const_1_minus_Cos_4_phi_data ) ) @@ -758,16 +781,19 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180)) + k0 + + +k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180)) + k4 * (1 - np.cos(4 * scanned_phi * np.pi / 180)) ) ) + elif use_const_1_plus_minus_Cos_x_values_bool is True: dihedral_energy = ( 1 / 2 * ( - +k3 * const_1_plus_Cos_3_phi_data + k0 + + +k3 * const_1_plus_Cos_3_phi_data + k4 * const_1_minus_Cos_4_phi_data ) ) @@ -778,7 +804,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k0 + + k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180)) ) ) @@ -788,7 +815,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * const_1_plus_Cos_1_phi_data + k0 * const_1_minus_Cos_0_phi_data + + k1 * const_1_plus_Cos_1_phi_data + k2 * const_1_minus_Cos_2_phi_data ) ) @@ -799,7 +827,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k0 + + +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180)) + k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180)) ) @@ -809,7 +838,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * const_1_plus_Cos_1_phi_data + k0 + + +k1 * const_1_plus_Cos_1_phi_data + k2 * const_1_minus_Cos_2_phi_data + k3 * const_1_plus_Cos_3_phi_data ) @@ -821,7 +851,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k0 + + +k1 * (1 + np.cos(1 * scanned_phi * np.pi / 180)) + k2 * (1 - np.cos(2 * scanned_phi * np.pi / 180)) + k3 * (1 + np.cos(3 * scanned_phi * np.pi / 180)) + k4 * (1 - np.cos(4 * scanned_phi * np.pi / 180)) @@ -833,7 +864,8 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): 1 / 2 * ( - +k1 * const_1_plus_Cos_1_phi_data + k0 + + +k1 * const_1_plus_Cos_1_phi_data + k2 * const_1_minus_Cos_2_phi_data + k3 * const_1_plus_Cos_3_phi_data + k4 * const_1_minus_Cos_4_phi_data @@ -844,8 +876,6 @@ def opls_dihedral(cos_powers_phi_and_constants_data, k0, k1, k2, k3, k4): # periodic dihedral function using - - def periodic_dihedral_n_1_2_3_4_5( phi_data, K_0, @@ -964,6 +994,78 @@ def RB_torsion_n_1_2_3_4_5(phi_data, k_0, k_1, k_2, k_3, k_4, k_5): return torsion_energy +# OPLS dihedral energy only function using +def opls_dihedral_n_1_2_3_4(phi_data, k_0, k_1, k_2, k_3, k_4): + """OPLS dihedral energy calculation with only the selected k-values. + + This is the OPLS dihedral energy calculation done with only the selected k-values. + However, all k-values are to be input as this makes entering the + k-values in a standard method for all the different OPLS dihderal + configurations, and more importantly allows it to be fit to the data + properly, which is not possible if all the k-values are included. + The k-values are in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...), + and the output ' dihedral_energy' energy are output in the same energy units. + + NOTE: ALL THE K-VALUE ENERGY UNITS MUST BE THE SAME. + + .. math:: + opls_dihedral_n_1 &= 1/2 *( + &= k0 + &= + k1 * (1 + cos(1 * phi)) + &= + k2 * (1 - cos(2 * phi)) + &= + k3 * (1 + cos(3 * phi)) + &= + k4 * (1 - cos(4 * phi)) + &= ) + + Parameters + ---------- + phi_data: float, int, or list of floats/int, in degrees + The 'phi_data' angle of the dihedral is in degrees + k_0: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) + The 'k0' value is the k-value for the opls dihedral where n=0, + or the constant without a cosine multiple. The 'k0' can be a + constant or zero (0). + user entered value, because it is not in the equation form. + k_1: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) + The 'k1' value is the k-value for the opls dihedral where n=1 + in the cosine multiple. + k_2: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) + The 'k2' value is the k-value for the opls dihedral where n=2 + in the cosine multiple. + NOTE: In this case, it is set to zero (0) regardless of the + user entered value, because it is not in the equation form. + k_3: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) + The 'k3' value is the k-value for the opls dihedral where n=3 + in the cosine multiple. + NOTE: In this case, it is set to zero (0) regardless of the + user entered value, because it is not in the equation form. + k_4: int, or float, in energy units (i.e., kcal/mol, kJ/mol, Kelvin, ...) + The 'k4' value is the k-value for the opls dihedral where n=4 + in the cosine multiple. + NOTE: In this case, it is set to zero (0) regardless of the + user entered value, because it is not in the equation form. + + Returns + ------- + dihedral_energy; float or list of floats, same as the k-value energy units + The OPLS dihedral energy from the valid phi_data and k-values + (i.e., k1 in this case). + """ + dihedral_energy = ( + 1 + / 2 + * ( + k_0 + + k_1 * (1 + np.cos(1 * (phi_data * np.pi / 180))) + + k_2 * (1 - np.cos(2 * (phi_data * np.pi / 180))) + + k_3 * (1 + np.cos(3 * (phi_data * np.pi / 180))) + + k_4 * (1 - np.cos(4 * (phi_data * np.pi / 180))) + ) + ) + + return dihedral_energy + + def get_r_squared(data_points, fitted_values): """Get the R**2 (R squared) values from the fitted equation.