diff --git a/cg_openmm/build/cg_build.py b/cg_openmm/build/cg_build.py index 890050d6..cc63f9c9 100755 --- a/cg_openmm/build/cg_build.py +++ b/cg_openmm/build/cg_build.py @@ -584,7 +584,7 @@ def check_force(cgmodel, force, force_type=None): openmm_nonbonded_energy.value_in_unit(unit.kilojoule_per_mole) ) * unit.kilojoule_per_mole - if energy_diff > 1E-1 * unit.kilojoule_per_mole: + if energy_diff > 1E0 * unit.kilojoule_per_mole: print("Error: The nonbonded potential energy computed by hand does not agree") print("with the value computed by OpenMM.") print(f"The value computed by OpenMM was: {openmm_nonbonded_energy}") diff --git a/cg_openmm/parameters/evaluate_energy.py b/cg_openmm/parameters/evaluate_energy.py index 2b4b929c..96e9f5af 100644 --- a/cg_openmm/parameters/evaluate_energy.py +++ b/cg_openmm/parameters/evaluate_energy.py @@ -746,7 +746,9 @@ def match_go_binary_interaction_energies(cgmodel, medoid_file, binary_interactio if type(binary_interaction_list) in [int, float]: # Single parameter - convert to list: - binary_interaction_list = list(binary_interaction_list) + binary_interaction_val = binary_interaction_list + binary_interaction_list = [] + binary_interaction_list.append(binary_interaction_val) for kappa in binary_interaction_list: # Create parameter update instruction dictionary @@ -782,8 +784,6 @@ def optimize_eps_scaling(x, param_dict, cgmodel, medoid_file, medoid_ref_energy, eps_scaled = particle["epsilon"]*x param_dict[f"{type_name}_epsilon"] = eps_scaled - print(f"param_dict: {param_dict}") - # Evaluate the energy with updated parameters: U_scaled, simulation = eval_energy(cgmodel, medoid_file, None, param_dict, verbose=verbose) @@ -795,6 +795,308 @@ def optimize_eps_scaling(x, param_dict, cgmodel, medoid_file, medoid_ref_energy, return U_diff +def optimize_go_sequences(cgmodel, medoid_file, native_contact_list, kappa=0.5, + n_restypes=2, seq_prepend=(0,)): + """ + Given a cgmodel with a topology and system, a reference medoid file, and native contact definitions, + determine an optimal sequence from a 'n_restypes' alphabet which minimizes the energy from native + contact pairs and maximizes the energy from non-native pairs. Energy of the native structure is kept constant + by scaling up all epsilon parameters to compensate for the use of a binary interaction parameter which scales down + the non-native interactions. For n_restypes=2, a brute force approach is used to enumerate over all possible sequences. + For larger n_restypes, SciPy differential evolution is used. There are no constraints on the relative amounts of each bead type. + + :param cgmodel: CGModel() class object to evaluate energy with. Note that backbone and sidechain bead types should begin with 'b', and 's', respectively. + :type cgmodel: class + + :param medoid_file: path to reference medoid file for full interaction energy target + :type medoid_file: str + + :param kappa: Binary interaction parameter for cross interactions (i.e. A-B residues), which are scaled by (1-kappa). + :type kappa: int + + :param n_restypes: sequence alphabet size (default=2) + :type n_restypes: int + + :param seq_prepend: If not an empty tuple, fix the specified sequence in place, and scan the rest of the residues (default=(0,)). Note that the first residue can always be fixed. + :type seq_prepend: tuple + + :returns: + - U_eval_diff_out - dictionary mapping sequences to (U_native-U_nonnative)/(U_native+U_nonnative) + - U_eval_native_out - dictionary mapping sequences to native energy + - U_eval_nonnative_out - dictionary mapping sequences to nonnative energy + """ + + # ***TODO: Check that the initial system is a homopolymer, and generalize the bead type parsing. + + # Get number of particles per monomer: + cgmodel_mono0 = cgmodel.monomer_types[0] + n_particle_per_mono = len(cgmodel_mono0["particle_sequence"]) + + # Get the number of beads: + num_beads = cgmodel.get_num_beads() + + # Get the number of monomers: + num_monomers = int(num_beads/n_particle_per_mono) + + # Construct the list of sequences: + seq_unique = [] + + prod_list = [] + for i in range(n_restypes): + prod_list.append(i) + + for seq in list(product(prod_list,repeat=(num_monomers-len(seq_prepend)))): + seq_full = seq_prepend+seq + seq_unique.append(seq_full) + + # These checks are very slow - usually it is faster to just compute all sequences: + + # if tuple(reversed(seq)) not in seq_unique: # No reverse duplicates + # We also want to avoid sequences that would be the same if 0 and 1 are swapped. + # For example 0010001 and 1101110 will give the same solution. + # For num_monomers > 2, there are more cases or repeated results, such as 001000 and 002000. + + # if n_restypes == 2: + # seq_array = np.array(seq) + # seq_swap = (seq_array*-1)+1 + # if tuple(seq_swap) not in seq_unique and tuple(reversed(seq_swap)) not in seq_unique: + # # These seq are tuples + # seq_unique.append(seq) + + # if n_restypes == 3: + # if 2 in seq: + # if 0 in seq and 1 in seq: + # # Adding a third type is only a new result if 0 and 1 are present. + # seq_unique.append(seq) + # else: + # seq_unique.append(seq) + # else: + # seq_array = np.array(seq) + # seq_swap = (seq_array*-1)+1 + # if tuple(seq_swap) not in seq_unique and tuple(reversed(seq_swap)) not in seq_unique: + # # These seq are tuples + # seq_unique.append(seq) + + print(f'scanning {len(seq_unique)} sequences') + + # The cgmodel should be set up as a Go model, with each particle its own type. + # These should be named '{type_name}{res_id}' + # This way, we implement the binary interaction parameters for each sequence. + + # Get full list of particle type names in order of bead index: + particle_type_name_all = [] + for bead_id in range(num_beads): + particle_type_name_all.append(cgmodel.get_particle_type_name(bead_id)) + + U_eval_diff_out = {} + U_eval_native_out = {} + U_eval_nonnative_out = {} + + # Get the nonbonded interaction list from cgmodel: + nonbonded_interaction_list = cgmodel.get_nonbonded_interaction_list() + nonbonded_exclusion_list = cgmodel.get_nonbonded_exclusion_list() + + nonnative_pair_list = [] + + for pair in nonbonded_interaction_list: + par1 = pair[0] + par2 = pair[1] + if ([par1,par2] not in nonbonded_exclusion_list and [par2,par1] not in nonbonded_exclusion_list and \ + [par1,par2] not in native_contact_list and [par2,par1] not in native_contact_list): + nonnative_pair_list.append(pair) + + # If different LJ parameters, we need to separate by pair type: + bb_bb_native_list = [] + bb_sc_native_list = [] + sc_sc_native_list = [] + + bb_bb_nonnative_list = [] + bb_sc_nonnative_list = [] + sc_sc_nonnative_list = [] + + # Also get the interaction types for each pair: + # The type names may be different for every particle. Just use the first letter + # 'b' for backbone and 's' for sidechain. + + for pair in native_contact_list: + type1 = particle_type_name_all[pair[0]][0:1] + type2 = particle_type_name_all[pair[1]][0:1] + if type1 == 'b' and type2 == 'b': + bb_bb_native_list.append(pair) + elif type1 == 's' and type2 == 's': + sc_sc_native_list.append(pair) + elif (type1 == 'b' and type2 == 's') or (type1 == 's' and type2 == 'b'): + bb_sc_native_list.append(pair) + + for pair in nonnative_pair_list: + type1 = particle_type_name_all[pair[0]][0:1] + type2 = particle_type_name_all[pair[1]][0:1] + if type1 == 'b' and type2 == 'b': + bb_bb_nonnative_list.append(pair) + elif type1 == 's' and type2 == 's': + sc_sc_nonnative_list.append(pair) + elif (type1 == 'b' and type2 == 's') or (type1 == 's' and type2 == 'b'): + bb_sc_nonnative_list.append(pair) + + # Load the native structure file: + if medoid_file[-3:] == 'dcd': + medoid_traj = md.load(medoid_file,top=md.Topology.from_openmm(cgmodel.topology)) + else: + medoid_traj = md.load(medoid_file) + + # Precompute the LJ functions (native and non-native) + # The arrays from mdtraj are [n_frames x n_pairs] + bb_bb_native_dist_array = md.compute_distances(medoid_traj,bb_bb_native_list) + bb_sc_native_dist_array = md.compute_distances(medoid_traj,bb_sc_native_list) + sc_sc_native_dist_array = md.compute_distances(medoid_traj,sc_sc_native_list) + + bb_bb_nonnative_dist_array = md.compute_distances(medoid_traj,bb_bb_nonnative_list) + bb_sc_nonnative_dist_array = md.compute_distances(medoid_traj,bb_sc_nonnative_list) + sc_sc_nonnative_dist_array = md.compute_distances(medoid_traj,sc_sc_nonnative_list) + + # Get the sigma and epsilon parameters: + bb_found = False + sc_found = False + for i in range(num_beads): + type_name_i = cgmodel.get_particle_type_name(i) + if type_name_i[0:1] == 'b': + bb_found = True + sigma_bb = cgmodel.get_particle_sigma(i).value_in_unit(unit.nanometer) + epsilon_bb = cgmodel.get_particle_epsilon(i).value_in_unit(unit.kilojoule_per_mole) + elif type_name_i[0:1] == 's': + sc_found = True + sigma_sc = cgmodel.get_particle_sigma(i).value_in_unit(unit.nanometer) + epsilon_sc = cgmodel.get_particle_epsilon(i).value_in_unit(unit.kilojoule_per_mole) + + if bb_found and sc_found: + break + + sigma_bb_sc = (sigma_bb+sigma_sc)/2 + epsilon_bb_sc = np.sqrt(epsilon_bb*epsilon_sc) + + LJ_12_factor_bb_bb = 4*epsilon_bb*np.power(sigma_bb,12) + LJ_6_factor_bb_bb = -4*epsilon_bb*np.power(sigma_bb,6) + + LJ_12_factor_bb_sc = 4*epsilon_bb_sc*np.power(sigma_bb_sc,12) + LJ_6_factor_bb_sc = -4*epsilon_bb_sc*np.power(sigma_bb_sc,6) + + LJ_12_factor_sc_sc = 4*epsilon_sc*np.power(sigma_sc,12) + LJ_6_factor_sc_sc = -4*epsilon_sc*np.power(sigma_sc,6) + + # These are the full LJ energies for all pair types split into native and non-native contributions, + # with no binary interaction scaling. Multiply element-wise by kappa vector specific to each sequence. + + LJ_full_bb_bb_native = LJ_12_factor_bb_bb/np.power(bb_bb_native_dist_array,12) + LJ_6_factor_bb_bb/np.power(bb_bb_native_dist_array,6) + LJ_full_bb_sc_native = LJ_12_factor_bb_sc/np.power(bb_sc_native_dist_array,12) + LJ_6_factor_bb_sc/np.power(bb_sc_native_dist_array,6) + LJ_full_sc_sc_native = LJ_12_factor_sc_sc/np.power(sc_sc_native_dist_array,12) + LJ_6_factor_sc_sc/np.power(sc_sc_native_dist_array,6) + + LJ_full_bb_bb_nonnative = LJ_12_factor_bb_bb/np.power(bb_bb_nonnative_dist_array,12) + LJ_6_factor_bb_bb/np.power(bb_bb_nonnative_dist_array,6) + LJ_full_bb_sc_nonnative = LJ_12_factor_bb_sc/np.power(bb_sc_nonnative_dist_array,12) + LJ_6_factor_bb_sc/np.power(bb_sc_nonnative_dist_array,6) + LJ_full_sc_sc_nonnative = LJ_12_factor_sc_sc/np.power(sc_sc_nonnative_dist_array,12) + LJ_6_factor_sc_sc/np.power(sc_sc_nonnative_dist_array,6) + + # For each sequence, multiply the LJ vectors by vectors containing 1 for same-type interaction, (1-kappa) for cross interactions (i.e., A-B) + + # Only need to initialize these once: + bb_bb_native_kappa_weights = np.zeros(len(bb_bb_native_list)) + bb_sc_native_kappa_weights = np.zeros(len(bb_sc_native_list)) + sc_sc_native_kappa_weights = np.zeros(len(sc_sc_native_list)) + + bb_bb_nonnative_kappa_weights = np.zeros(len(bb_bb_nonnative_list)) + bb_sc_nonnative_kappa_weights = np.zeros(len(bb_sc_nonnative_list)) + sc_sc_nonnative_kappa_weights = np.zeros(len(sc_sc_nonnative_list)) + + n = 0 + for seq in seq_unique: + + # Sequence by residue: + seq_array_res = np.asarray(seq) + + # Sequence by particle: + seq_array_par = np.zeros(num_beads) + i = 0 + for res_seq in seq_array_res: + for par in range(n_particle_per_mono): + seq_array_par[i] = res_seq + i += 1 + + # Construct binary interaction weighting vectors: + # Unlike residue interactions get (1-kappa) weight, like residues get 1 weight. + + # native bb-bb + i = 0 + for pair in bb_bb_native_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + bb_bb_native_kappa_weights[i] = 1 + else: + bb_bb_native_kappa_weights[i] = (1-kappa) + i += 1 + + # native bb-sc + i = 0 + for pair in bb_sc_native_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + bb_sc_native_kappa_weights[i] = 1 + else: + bb_sc_native_kappa_weights[i] = (1-kappa) + i += 1 + + # native sc-sc + i = 0 + for pair in sc_sc_native_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + sc_sc_native_kappa_weights[i] = 1 + else: + sc_sc_native_kappa_weights[i] = (1-kappa) + i += 1 + + # nonnative bb-bb + i = 0 + for pair in bb_bb_nonnative_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + bb_bb_nonnative_kappa_weights[i] = 1 + else: + bb_bb_nonnative_kappa_weights[i] = (1-kappa) + i += 1 + + # nonnative bb-sc + i = 0 + for pair in bb_sc_nonnative_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + bb_sc_nonnative_kappa_weights[i] = 1 + else: + bb_sc_nonnative_kappa_weights[i] = (1-kappa) + i += 1 + + # nonnative sc-sc + i = 0 + for pair in sc_sc_nonnative_list: + if seq_array_par[pair[0]] == seq_array_par[pair[1]]: + sc_sc_nonnative_kappa_weights[i] = 1 + else: + sc_sc_nonnative_kappa_weights[i] = (1-kappa) + i += 1 + + # Compute the native LJ energy of current sequence: + native_energy_total = ( + (LJ_full_bb_bb_native*bb_bb_native_kappa_weights).sum() + \ + (LJ_full_bb_sc_native*bb_sc_native_kappa_weights).sum() + \ + (LJ_full_sc_sc_native*sc_sc_native_kappa_weights).sum() + ) + + # Compute the nonnative LJ energy of current sequence: + nonnative_energy_total = ( + (LJ_full_bb_bb_nonnative*bb_bb_nonnative_kappa_weights).sum() + \ + (LJ_full_bb_sc_nonnative*bb_sc_nonnative_kappa_weights).sum() + \ + (LJ_full_sc_sc_nonnative*sc_sc_nonnative_kappa_weights).sum() + ) + + U_eval_diff_out[f'{seq}'] = (native_energy_total - nonnative_energy_total) / (native_energy_total + nonnative_energy_total) + U_eval_native_out[f'{seq}'] = native_energy_total + U_eval_nonnative_out[f'{seq}'] = nonnative_energy_total + + return U_eval_diff_out, U_eval_native_out, U_eval_nonnative_out + + def eval_energy_sequences(cgmodel, file_list, temperature_list, monomer_list, sequence=None, output_data='output/output.nc', num_intermediate_states=3, n_trial_boot=200, plot_dir='', frame_begin=0, frame_end=-1, sample_spacing=1, sparsify_stride=1, verbose=False, n_cpu=1): diff --git a/cg_openmm/tests/test_data/rebuild_test_cgmodels.py b/cg_openmm/tests/test_data/rebuild_test_cgmodels.py index 43a65a51..e7b442f1 100644 --- a/cg_openmm/tests/test_data/rebuild_test_cgmodels.py +++ b/cg_openmm/tests/test_data/rebuild_test_cgmodels.py @@ -5,13 +5,16 @@ import numpy as np import openmm from cg_openmm.cg_model.cgmodel import CGModel +import mdtraj as md from openmm import unit from openmm.app.pdbfile import PDBFile +from cg_openmm.parameters.secondary_structure import get_native_contacts output_file1 = "stored_cgmodel.pkl" output_file2 = "stored_cgmodel_per1_3.pkl" output_file3 = "linear_24mer/stored_cgmodel_24mer_linear.pkl" output_file4 = "stored_cgmodel_binary_interaction.pkl" +output_file5 = "stored_cgmodel_helix_3sc_open_triangle_8mer.pkl" def get_test_cmgodel_1_1_helix(output_file): # Generate a cgmodel with 1 backbone bead, 1 sidechain bead per residue, @@ -422,7 +425,6 @@ def get_test_cmgodel_1_1_helix_binary_interaction(output_file): # Residue sequence: sequence = 24 * [A] - #-------------------------------# # Binary interaction parameters # #-------------------------------# @@ -505,9 +507,239 @@ def get_test_cmgodel_1_1_helix_binary_interaction(output_file): return cgmodel + +def get_test_cmgodel_3sc_open_triangle_Go(output_file): + # Generate a cgmodel with 1 backbone bead, 3 sidechain bead per residue, + # in open triangle configuration. + + n_mono = 8 + + include_bond_forces = True + include_bond_angle_forces = True + include_nonbonded_forces = True + include_torsion_forces = True + constrain_bonds = False + + # Set backbone-backbone exclusions to [0,0,1], all others to [1,1,1] + exclusions = { + "default_exclusions": [1,1,1], + } + + for i in range(n_mono): + for j in range(i+1,n_mono): + exclusions[f"b{i}_b{j}_exclusions"] = [0,0,1] + + #--------------------------------------------# + # Particle definitions and oligomer topology # + #--------------------------------------------# + + mass = 100.0 * unit.amu + + # All residues get the same sigma, epsilon but different particle type names + particle_type_list = [] + for i in range(n_mono): + particle_dict_b = { + "particle_type_name": f"b{i}", + "sigma": 2.25 * unit.angstrom, + "epsilon": 1.5 * unit.kilojoules_per_mole, + "mass": mass + } + + particle_dict_sa = { + "particle_type_name": f"sa{i}", + "sigma": 2.007042724 * unit.angstrom, + "epsilon": 5.0 * unit.kilojoules_per_mole, + "mass": mass + } + + particle_dict_sb = { + "particle_type_name": f"sb{i}", + "sigma": 2.007042724 * unit.angstrom, + "epsilon": 5.0 * unit.kilojoules_per_mole, + "mass": mass + } + + particle_dict_sc = { + "particle_type_name": f"sc{i}", + "sigma": 2.007042724 * unit.angstrom, + "epsilon": 5.0 * unit.kilojoules_per_mole, + "mass": mass + } + + exec(f'b{i} = particle_dict_b') + exec(f'sa{i} = particle_dict_sa') + exec(f'sb{i} = particle_dict_sb') + exec(f'sc{i} = particle_dict_sc') + + particle_type_list.append(eval(f'b{i}')) + particle_type_list.append(eval(f'sa{i}')) + particle_type_list.append(eval(f'sb{i}')) + particle_type_list.append(eval(f'sc{i}')) + + # Monomer definition + sequence = [] + monomer_types = [] + + for i in range(n_mono): + monomer_dict = { + "monomer_name": f"A{i}", + "particle_sequence": [eval(f"b{i}"),eval(f"sa{i}"),eval(f"sb{i}"),eval(f"sc{i}")], + "bond_list": [[0,1],[1,2],[2,3]], # open-form triangle + "start": 0, + "end": 0, + } + + exec(f'A{i} = monomer_dict') + + sequence.append(eval(f'A{i}')) + monomer_types.append(eval(f'A{i}')) + + #--------------------------# + # Harmonic bond parameters # + #--------------------------# + + # Bond definitions + # Set the bb-bb types to be the default, and loop over the bb-sc types: + bond_lengths = {"default_bond_length": 0.250660083 * unit.nanometer} + + for i in range(n_mono): + bond_lengths[f"b{i}_sa{i}_bond_length"] = 0.238918445 * unit.nanometer + bond_lengths[f"sa{i}_sb{i}_bond_length"] = 0.225282929 * unit.nanometer + bond_lengths[f"sa{i}_sc{i}_bond_length"] = 0.225282929 * unit.nanometer + bond_lengths[f"sb{i}_sc{i}_bond_length"] = 0.225282929 * unit.nanometer + + bond_force_constants = { + "default_bond_force_constant": 15000 * unit.kilojoule_per_mole / unit.nanometer / unit.nanometer + } + + #---------------------------# + # Harmonic angle parameters # + #---------------------------# + + # Bond angle definitions + bond_angle_force_constants = { + "default_bond_angle_force_constant": 0 * unit.kilojoule_per_mole / unit.radian / unit.radian, + } + for i in range(n_mono-2): + bond_angle_force_constants[f"b{i}_b{i+1}_b{i+2}_bond_angle_force_constant"] = 150 * unit.kilojoule_per_mole / unit.radian / unit.radian + + # Set the bb-bb-sc bond angles to be the default, loop over the bb-bb-bb types: + equil_bond_angles = {} + for i in range(n_mono-2): + equil_bond_angles[f"b{i}_b{i+1}_b{i+2}_equil_bond_angle"] = 83.03316 * unit.degrees + + #-----------------------------# + # Periodic torsion parameters # + #-----------------------------# + + # Torsion angle definitions + torsion_force_constants = {"default_torsion_force_constant": 0.0 * unit.kilojoule_per_mole} + for i in range(n_mono-3): + torsion_force_constants[f"b{i}_b{i+1}_b{i+2}_b{i+3}_torsion_force_constant"] = 0 * unit.kilojoule_per_mole + + # Only backbone types get non-zero force constants, so can set a single default for these: + torsion_phase_angles = { + "default_torsion_phase_angle": (33.434975-180) * unit.degrees, + } + + torsion_periodicities = { + "default_torsion_periodicity": 1, + } + + #------------------------------------------# + # Binary interaction parameters (Go model) # + #------------------------------------------# + + native_positions_file = "../test_structures/helix_3sc_open_triangle_8mer.pdb" + positions = PDBFile(native_positions_file).getPositions() + + # First, build a dummy cgmodel with no binary interaction parameters: + cgmodel_non_go = CGModel( + particle_type_list=particle_type_list, + bond_lengths=bond_lengths, + bond_force_constants=bond_force_constants, + bond_angle_force_constants=bond_angle_force_constants, + torsion_force_constants=torsion_force_constants, + equil_bond_angles=equil_bond_angles, + torsion_phase_angles=torsion_phase_angles, + torsion_periodicities=torsion_periodicities, + include_nonbonded_forces=include_nonbonded_forces, + include_bond_forces=include_bond_forces, + include_bond_angle_forces=include_bond_angle_forces, + include_torsion_forces=include_torsion_forces, + constrain_bonds=constrain_bonds, + exclusions=exclusions, + positions=positions, + sequence=sequence, + monomer_types=monomer_types, + ) + + native_contact_distance_cutoff = 3.15 * unit.angstrom + + # Determine all native contacts: + native_contact_list, native_contact_distances, contact_type_dict = get_native_contacts( + cgmodel_non_go, + native_positions_file, + native_contact_distance_cutoff, + ) + + # Beads per residue: + n_per_res = 4 + + # Binary interaction parameters: + # (any not set are assigned kappa=0) + binary_interaction_parameters = {} + + # Get the particle type name lists: + type_name_list = [] + for i in range(n_mono): + type_name_list.append('b') + type_name_list.append('sa') + type_name_list.append('sb') + type_name_list.append('sc') + + # The beads get ordered [b0,s0,b1,s1,...,b23,s23] + for i in range(n_mono*n_per_res): + for j in range(i+1,n_mono*n_per_res): + if ([i,j] not in native_contact_list and [j,i] not in native_contact_list): + # If not in native contact list, scale down the interaction by (1-kappa_go) + + type_i = type_name_list[i] + type_j = type_name_list[j] + + resid_i = i//n_per_res + resid_j = j//n_per_res + + binary_interaction_parameters[f"{type_i}{resid_i}_{type_j}{resid_j}_binary_interaction"] = 0.99 + + # Build a coarse grained model + cgmodel = CGModel( + particle_type_list=particle_type_list, + bond_lengths=bond_lengths, + bond_force_constants=bond_force_constants, + bond_angle_force_constants=bond_angle_force_constants, + torsion_force_constants=torsion_force_constants, + equil_bond_angles=equil_bond_angles, + torsion_phase_angles=torsion_phase_angles, + torsion_periodicities=torsion_periodicities, + binary_interaction_parameters=binary_interaction_parameters, + include_nonbonded_forces=include_nonbonded_forces, + include_bond_forces=include_bond_forces, + include_bond_angle_forces=include_bond_angle_forces, + include_torsion_forces=include_torsion_forces, + constrain_bonds=constrain_bonds, + exclusions=exclusions, + positions=positions, + sequence=sequence, + monomer_types=monomer_types, + ) + + # store the cg model so that we can do various analyses. + cgmodel.export(output_file) # Generate the test cgmodels: cgmodel1 = get_test_cmgodel_1_1_helix(output_file1) cgmodel2 = get_test_cmgodel_1_1_helix_per1_3(output_file2) cgmodel3 = get_test_cmgodel_linear(output_file3) -cgmodel4 = get_test_cmgodel_1_1_helix_binary_interaction(output_file4) \ No newline at end of file +cgmodel4 = get_test_cmgodel_1_1_helix_binary_interaction(output_file4) +cgmodel5 = get_test_cmgodel_3sc_open_triangle_Go(output_file5) \ No newline at end of file diff --git a/cg_openmm/tests/test_data/stored_cgmodel_helix_3sc_open_triangle_8mer.pkl b/cg_openmm/tests/test_data/stored_cgmodel_helix_3sc_open_triangle_8mer.pkl new file mode 100644 index 00000000..47127301 Binary files /dev/null and b/cg_openmm/tests/test_data/stored_cgmodel_helix_3sc_open_triangle_8mer.pkl differ diff --git a/cg_openmm/tests/test_energy_eval.py b/cg_openmm/tests/test_energy_eval.py index f4eca087..4241c932 100644 --- a/cg_openmm/tests/test_energy_eval.py +++ b/cg_openmm/tests/test_energy_eval.py @@ -12,6 +12,7 @@ from cg_openmm.parameters.evaluate_energy import * from cg_openmm.parameters.reweight import (get_opt_temperature_list, get_temperature_list) +from cg_openmm.parameters.secondary_structure import get_native_contacts from cg_openmm.thermo.calc import * from numpy.testing import assert_allclose, assert_almost_equal from openmm import unit @@ -1204,12 +1205,43 @@ def test_optimize_epsilon_scaling_go_models(tmpdir): cgmodel, medoid_file, binary_interaction_list, - verbose=True, + verbose=False, ) # Check that all scaling factors are greater than 1: for kappa in binary_interaction_list: - assert eps_scaling['kappa_{kappa}'] >= 1.0 + assert eps_scaling[f'kappa_{kappa}'] >= 1.0 + + +def test_native_sequence_optimization_AB_pdb(tmpdir): + """ + Test sequence optimization for 3sc triangle model, where the target is minimizing the + contribution of non-native contacts to the energy of the native folded structure. + (native structure is pdb file) + """ + + # Path to native structure + native_structure_file = f"{structures_path}/helix_3sc_open_triangle_8mer.pdb" + + # Load in cgmodel + cgmodel = pickle.load(open(f"{data_path}/stored_cgmodel_helix_3sc_open_triangle_8mer.pkl", "rb" )) + + # Get native contact list: + native_contact_distance_cutoff = 3.15 * unit.angstrom + native_contact_list, native_contact_distances, contact_type_dict = get_native_contacts( + cgmodel, + native_structure_file, + native_contact_distance_cutoff, + ) + + U_eval_diff_out, U_eval_native_out, U_eval_nonnative_out = optimize_go_sequences( + cgmodel, native_structure_file, native_contact_list, kappa=0.99, n_restypes=2, + seq_prepend=(0,), + ) + + assert len(U_eval_diff_out) == 2**7 + assert len(U_eval_native_out) == 2**7 + assert len(U_eval_nonnative_out) == 2**7 def test_reeval_heat_capacity(tmpdir): diff --git a/cg_openmm/tests/test_structures/helix_3sc_open_triangle_8mer.pdb b/cg_openmm/tests/test_structures/helix_3sc_open_triangle_8mer.pdb new file mode 100644 index 00000000..44b7b33a --- /dev/null +++ b/cg_openmm/tests/test_structures/helix_3sc_open_triangle_8mer.pdb @@ -0,0 +1,65 @@ +ATOM 1 b01 A0 A 1 1.572 0.000 0.000 1.00 0.00 +ATOM 2 sa02 A0 A 1 3.576 -1.225 0.440 1.00 0.00 +ATOM 3 sb03 A0 A 1 3.576 0.993 0.839 1.00 0.00 +ATOM 4 sc04 A0 A 1 3.576 0.231 -1.280 1.00 0.00 +ATOM 5 b11 A1 A 2 -0.305 1.542 0.619 1.00 0.00 +ATOM 6 sa12 A1 A 2 0.506 3.746 1.059 1.00 0.00 +ATOM 7 sb13 A1 A 2 -1.668 3.315 1.459 1.00 0.00 +ATOM 8 sc14 A1 A 2 -0.920 3.463 -0.662 1.00 0.00 +ATOM 9 b21 A2 A 3 -1.454 -0.599 1.237 1.00 0.00 +ATOM 10 sa22 A2 A 3 -3.773 -0.229 1.677 1.00 0.00 +ATOM 11 sb23 A2 A 3 -2.928 -2.279 2.076 1.00 0.00 +ATOM 12 sc24 A2 A 3 -3.219 -1.573 -0.042 1.00 0.00 +ATOM 13 b31 A3 A 4 0.869 -1.309 1.856 1.00 0.00 +ATOM 14 sa32 A3 A 4 0.957 -3.656 2.297 1.00 0.00 +ATOM 15 sb33 A3 A 4 2.804 -2.430 2.696 1.00 0.00 +ATOM 16 sc34 A3 A 4 2.169 -2.853 0.575 1.00 0.00 +ATOM 17 b41 A4 A 5 1.116 1.107 2.474 1.00 0.00 +ATOM 18 sa42 A4 A 5 3.402 1.649 2.914 1.00 0.00 +ATOM 19 sb43 A4 A 5 1.841 3.222 3.314 1.00 0.00 +ATOM 20 sc44 A4 A 5 2.378 2.682 1.195 1.00 0.00 +ATOM 21 b51 A5 A 6 -1.302 0.880 3.093 1.00 0.00 +ATOM 22 sa52 A5 A 6 -2.276 3.016 3.534 1.00 0.00 +ATOM 23 sb53 A5 A 6 -3.519 1.180 3.933 1.00 0.00 +ATOM 24 sc54 A5 A 6 -3.092 1.812 1.812 1.00 0.00 +ATOM 25 b61 A6 A 7 -0.611 -1.448 3.711 1.00 0.00 +ATOM 26 sa62 A6 A 7 -2.518 -2.818 4.151 1.00 0.00 +ATOM 27 sb63 A6 A 7 -0.475 -3.681 4.551 1.00 0.00 +ATOM 28 sc64 A6 A 7 -1.179 -3.384 2.432 1.00 0.00 +ATOM 29 b71 A7 A 8 1.539 -0.319 4.331 1.00 0.00 +ATOM 30 sa72 A7 A 8 3.254 -1.924 4.771 1.00 0.00 +ATOM 31 sb73 A7 A 8 3.703 0.248 5.170 1.00 0.00 +ATOM 32 sc74 A7 A 8 3.548 -0.500 3.049 1.00 0.00 +TER +CONECT 1 2 +CONECT 2 3 +CONECT 3 4 +CONECT 1 5 +CONECT 5 6 +CONECT 6 7 +CONECT 7 8 +CONECT 5 9 +CONECT 9 10 +CONECT 10 11 +CONECT 11 12 +CONECT 9 13 +CONECT 13 14 +CONECT 14 15 +CONECT 15 16 +CONECT 13 17 +CONECT 17 18 +CONECT 18 19 +CONECT 19 20 +CONECT 17 21 +CONECT 21 22 +CONECT 22 23 +CONECT 23 24 +CONECT 21 25 +CONECT 25 26 +CONECT 26 27 +CONECT 27 28 +CONECT 25 29 +CONECT 29 30 +CONECT 30 31 +CONECT 31 32 +END diff --git a/cg_openmm/utilities/helix_optimize_geometry.py b/cg_openmm/utilities/helix_optimize_geometry.py index 844fc1b5..42b969d4 100644 --- a/cg_openmm/utilities/helix_optimize_geometry.py +++ b/cg_openmm/utilities/helix_optimize_geometry.py @@ -70,7 +70,7 @@ def optimize_helix_simple(n_particle_bb, sigma, epsilon, sidechain=True, DE_pops r_eq = sigma*np.power(2,(1/6)) # Get particle positions: - xyz_par = get_helix_coordinates(r_opt,c_opt,t_par) + xyz_par = get_helix_backbone_coordinates(r_opt,c_opt,t_par) if sidechain: # Place sidechain particles normal to helix with same bond length as bb_bb @@ -288,7 +288,7 @@ def optimize_helix_openmm_energy(n_particle_bb, sigma_bb, sigma_sc, epsilon_bb, r_eq_sc = sigma_sc_val*np.power(2,(1/6)) # Get particle positions: - xyz_par = get_helix_coordinates(r_opt,c_opt,t_par) + xyz_par = get_helix_backbone_coordinates(r_opt,c_opt,t_par) # Place sidechain particles normal to helix if bond_dist_sc == None: @@ -382,7 +382,7 @@ def compute_LJ_helix_energy_simple(geo, sigma, epsilon, n_particle_bb, sidechain for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # Add any sidechain beads if sidechain: @@ -428,7 +428,7 @@ def compute_LJ_helix_openmm_energy(geo, simulation, bb_array, sc_array, n_partic for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -478,7 +478,7 @@ def compute_LJ_helix_openmm_energy_constrained( for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just diff --git a/cg_openmm/utilities/helix_optimize_nonbonded.py b/cg_openmm/utilities/helix_optimize_nonbonded.py index 8a96ac93..ea85e75b 100644 --- a/cg_openmm/utilities/helix_optimize_nonbonded.py +++ b/cg_openmm/utilities/helix_optimize_nonbonded.py @@ -8,11 +8,11 @@ from openmm import LangevinIntegrator, unit from openmm.app import Simulation from openmm.app.pdbfile import PDBFile -from scipy.optimize import differential_evolution, root_scalar, brute +from scipy.optimize import differential_evolution, root_scalar, brute, LinearConstraint -def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, - bond_dist_bb=None, bond_dist_sc=None, equal_bonds=True, DE_popsize=50, +def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, epsilon_sc=None, + bond_dist_bb=None, bond_dist_sc=None, equal_bonds=True, sigma_bb=None, DE_popsize=50, pdbfile='LJ_helix_openmm_energy.pdb', plotfile='LJ_helix_openmm_energy.pdf', exclusions={}): """ 1sc model: Optimize backbone and sidechain particle parameters along a helix with specified radius and @@ -27,14 +27,17 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, :param n_particle_bb: Number of backbone particles to model :type n_particle_bb: int - :param bond_dist_bb: bond distance for bb-bb bonds. If None, bond distance will also be optimized. + :param bond_dist_bb: bond distance for bb-bb bonds. If None, bond distance will also be optimized. (default=None) :type bond_dist_bb: Quantity - :param bond_dist_sc: bond distance for bb-sc bonds. If None, bond distance will also be optimized. + :param bond_dist_sc: bond distance for bb-sc bonds. If None, bond distance will also be optimized. If 'LJ', the equilibrium nonbonded distance will be used instead. (default=None). :type bond_dist_sc: Quantity :param equal_bonds: option to constrain bb-sc bond distance to equal bb-bb bond distance. If True, any specified bond distances are ignored. (default=True) :type equal_bonds: bool + + :param sigma_bb: sigma LJ parameter for backbone. If None, sigma_bb will also be optimized. (default=None) + :type sigma_bb: Quantity :param DE_popsize: population size to use in SciPy differential_evolution solver (default=50) :type DE_popsize: int @@ -67,13 +70,14 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, # Set initial epsilon parameters epsilon_bb = 1.0 * unit.kilojoule_per_mole - epsilon_sc = 1.0 * unit.kilojoule_per_mole + if epsilon_sc is None: + epsilon_sc = 1.0 * unit.kilojoule_per_mole # Set initial sigma parameters - sigma_bb = 1.0 * unit.angstrom - sigma_sc = 1.0 * unit.angstrom + sigma_bb_init = 1.0 * unit.angstrom + sigma_sc_init = 1.0 * unit.angstrom - cgmodel = get_helix_cgmodel(sigma_bb,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) + cgmodel = get_helix_cgmodel(sigma_bb_init,sigma_sc_init,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) # Get particle type lists and bonded lists: (particle_type_list, bb_array, sc_array, bb_bond_list, sc_bond_list, @@ -86,9 +90,9 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, integrator = LangevinIntegrator( 0.0 * unit.kelvin, friction, simulation_time_step.in_units_of(unit.picosecond) ) - simulation = Simulation(cgmodel.topology, cgmodel.system, integrator) + simulation = Simulation(cgmodel.topology, cgmodel.system, integrator) - if bond_dist_bb is None and bond_dist_sc is None: + if (bond_dist_bb is None and bond_dist_sc is None) or equal_bonds: #-----------------------------# # Unconstrained bonds version # @@ -102,6 +106,8 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, if equal_bonds: # Set optimization bounds [t, sigma_bb, sigma_sc]: # Use a minimium of 3 residues/turn + + # TODO: add case where we specific sigma_bb and not bond lengths bounds = [(0.01,2*np.pi/3),(r/50,15*r),(r/50,15*r)] opt_sol = differential_evolution( @@ -134,18 +140,82 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, sigma_sc_opt = opt_sol.x[2] r_bs_opt = opt_sol.x[3] - else: + elif bond_dist_bb is not None and bond_dist_sc is None: - #---------------------------# - # Constrained bonds version # - #---------------------------# + #----------------------------# + # 1 constrained bond version # + #----------------------------# + + bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) + + params = (simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc, sigma_bb) + + if sigma_bb is None: + # Set optimization bounds [sigma_bb, sigma_sc, r_bs]: + bounds = [(r/50,15*r),(r/50,15*r),(r/50,15*r)] + else: + # Set optimization bounds [sigma_sc, r_bs]: + bounds = [(r/50,15*r),(r/50,15*r)] + + opt_sol = differential_evolution( + compute_helix_openmm_energy_vary_LJ_constrained, + bounds, + args=params, + polish=True, + popsize=DE_popsize, + ) + + if sigma_bb is None: + sigma_bb_opt = opt_sol.x[0] + sigma_sc_opt = opt_sol.x[1] + r_bs_opt = opt_sol.x[2] + else: + sigma_bb_opt = sigma_bb.value_in_unit(unit.angstrom) + sigma_sc_opt = opt_sol.x[0] + r_bs_opt = opt_sol.x[1] + + # Compute particle spacing based on bond constraints + t_delta_opt = get_t_from_bond_distance(r,c,bond_dist_bb) + if t_delta_opt < 0: + print(t_delta_opt) + t_delta_opt *= -1 + + + elif bond_dist_bb is None and bond_dist_sc is not None and bond_dist_sc != 'LJ': + + #----------------------------# + # 1 constrained bond version # + #----------------------------# + + bond_dist_sc = bond_dist_sc.value_in_unit(unit.angstrom) + + params = (simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc, sigma_bb) + + # Set optimization bounds [sigma_bb, sigma_sc, t]: + bounds = [(r/50,15*r),(r/50,15*r),(r/50,15*r)] + + opt_sol = differential_evolution( + compute_helix_openmm_energy_vary_LJ_constrained, + bounds, + args=params, + polish=True, + popsize=DE_popsize, + ) + + sigma_bb_opt = opt_sol.x[0] + sigma_sc_opt = opt_sol.x[1] + t_delta_opt = opt_sol.x[2] + + elif bond_dist_bb is not None and bond_dist_sc is not None and bond_dist_sc != 'LJ': - # For now, we have to specify both bb-bb and bb-sc bond distances + #-----------------------------# + # 2 constrained bonds version # + #-----------------------------# bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) bond_dist_sc = bond_dist_sc.value_in_unit(unit.angstrom) - params = (simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc) + params = (simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc, sigma_bb) # Set optimization bounds [sigma_bb, sigma_sc]: bounds = [(r/50,15*r),(r/50,15*r)] @@ -166,7 +236,48 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, if t_delta_opt < 0: print(t_delta_opt) t_delta_opt *= -1 + + elif bond_dist_bb is not None and bond_dist_sc == 'LJ': + # TODO: add the last case with bond_dist_bb is None and bond_dist_sc == 'LJ' + + #----------------------------# + # LJ distance bb-sc bonds # + #----------------------------# + + bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) + + params = (simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc, sigma_bb) + + if sigma_bb is None: + # Set optimization bounds [sigma_bb, sigma_sc]: + bounds = [(r/50,15*r),(r/50,15*r)] + else: + # Set optimization bounds [sigma_sc]: + bounds = [(r/50,15*r)] + opt_sol = differential_evolution( + compute_helix_openmm_energy_vary_LJ_constrained, + bounds, + args=params, + polish=True, + popsize=DE_popsize, + ) + + if sigma_bb is None: + sigma_bb_opt = opt_sol.x[0] + sigma_sc_opt = opt_sol.x[1] + r_bs_opt = (sigma_bb_opt+sigma_sc_opt)/2*np.power(2,(1/6)) + else: + sigma_bb_opt = sigma_bb.value_in_unit(unit.angstrom) + sigma_sc_opt = opt_sol.x[0] + r_bs_opt = r_bs_opt = (sigma_bb_opt+sigma_sc_opt)/2*np.power(2,(1/6)) + + # Compute particle spacing based on bond constraints + t_delta_opt = get_t_from_bond_distance(r,c,bond_dist_bb) + if t_delta_opt < 0: + print(t_delta_opt) + t_delta_opt *= -1 + # Determine backbone particle parametric coordinates: t_par = np.zeros(n_particle_bb) for i in range(n_particle_bb): @@ -177,7 +288,7 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, r_eq_sc = sigma_sc_opt*np.power(2,(1/6)) # Get particle positions: - xyz = get_helix_coordinates(r,c,t_par) + xyz = get_helix_backbone_coordinates(r,c,t_par) # Place sidechain particles normal to helix r_bb = dist_unitless(xyz[0,:],xyz[1,:]) @@ -185,6 +296,8 @@ def optimize_helix_LJ_parameters(radius, pitch, n_particle_bb, if equal_bonds: # Use optimized bond length from first two backbone beads: r_bs = r_bb + elif bond_dist_sc == 'LJ': + r_bs = r_bs_opt else: if bond_dist_sc is not None: # Use specified bb-sc bond distance: @@ -425,7 +538,7 @@ def optimize_helix_LJ_parameters_energy_diff(radius, pitch, n_particle_bb, t_par[i] = i*t_delta_opt # Get particle positions: - xyz = get_helix_coordinates(r,c,t_par) + xyz = get_helix_backbone_coordinates(r,c,t_par) # Place sidechain particles normal to helix r_bb = dist_unitless(xyz[0,:],xyz[1,:]) @@ -662,7 +775,7 @@ def optimize_helix_LJ_parameters_2sc(radius, pitch, n_particle_bb, sigma_bb, r_eq_sc2 = sigma_sc2_opt*np.power(2,(1/6)) # Get particle positions: - xyz = get_helix_coordinates(r,c,t_par) + xyz = get_helix_backbone_coordinates(r,c,t_par) # Place sidechain particles normal to helix: side_xyz1 = np.zeros((n_particle_bb,3)) @@ -755,9 +868,436 @@ def optimize_helix_LJ_parameters_2sc(radius, pitch, n_particle_bb, sigma_bb, return opt_sol, geometry +def optimize_helix_LJ_parameters_2sc_rotation(radius, pitch, n_particle_bb, sigma_bb, + bond_dist_bb, equal_sc=True, DE_popsize=200, pdbfile='LJ_helix_2sc_opt_rotation.pdb', + n_rotation_angles=2, alignment='center', exclusions={}): + """ + With specified radius, pitch, backbone sigma, and backbone-backbone bond length, + optimize the sidechain sigma, bb-sc, sc-sc bond lengths, and rotation angle of the + sidechain particle about the + + :param radius: fixed helical radius + :type radius: Quantity + + :param pitch: fixed helical pitch (c*2*pi) + :type pitch: Quantity + + :param n_particle_bb: Number of backbone particles to model + :type n_particle_bb: int + + :param sigma_bb: LJ sigma parameter for backbone beads. If None, sigma_bb and backbone particle spacing will also be optimized. + :type sigma_bb: Quantity + + :param bond_dist_bb: bond distance for bb-bb bonds. If None, sigma_bb and backbone particle spacing will also be optimized. + :type bond_dist_bb: Quantity + + :param equal_sc: Option for keeping the sigma_sc of each sidechain bead the same. If False, each sigma_sc will be varied independently. (default=True) + :type equal_sc: bool + + :param DE_popsize: population size to use in SciPy differential_evolution solver (default=50) + :type DE_popsize: int + + :param pdbfile: Path to pdb file for saving the helical structure (default='LJ_helix_2sc_opt.pdb') + :type pdbfile: str + + :param n_rotation_angles: number of independent sidechain rotation angles. For example 2 treats odd and even residue rotations independently. + :type n_rotation_angles: str + + :param alignment: sidechain alignment scheme - can be 'center' (center of sidechain group is fixed normal to backbone) or 'first' (first bead is normal to backbone) + :type alignment: str + + :param exclusions: pass cg_openmm exclusion rules to the cgmodel (by default [0,0,1] is applied to bb-bb, [0,1,1] to bb-sc, sc-sc) + :type exclusions: dict + + :returns: + - opt_sol - Results from scipy.optimize (dict) + - geometry - Dictionary containing key geometric parameters of the optimized helix + """ + + # Check input: + if alignment not in ['center','first']: + print(f'Error: invalid alignment parameter {alignment}') + exit() + + if alignment == 'center' and not equal_sc: + print(f'Error: center sidechain group alignment not compatible with non-equal sigma_sc') + exit() + + if n_rotation_angles < 1: + print(f'n_rotation_angles must be at least 1') + exit() + + r_unit = radius.unit + # Use angstrom for writing pdb file: + radius = radius.value_in_unit(unit.angstrom) + pitch = pitch.value_in_unit(unit.angstrom) + + r = radius + c = pitch/(2*np.pi) # Helical rise parameter + + # Here we need to create a cgmodel + + # Set initial epsilon parameters + epsilon_bb = 1.0 * unit.kilojoule_per_mole + epsilon_sc = 1.0 * unit.kilojoule_per_mole + + # Set initial sigma parameters + if sigma_bb is None: + sigma_bb_init = 1.0 * unit.angstrom + else: + sigma_bb_init = sigma_bb + sigma_sc = 1.0 * unit.angstrom + + if bond_dist_bb is None and sigma_bb is not None: + print(f'Warning: setting sigma_bb to None for variable backbone optimization') + sigma_bb = None + + if bond_dist_bb is not None and sigma_bb is None: + print(f'Warning: setting bond_dist_bb to None for variable backbone optimization') + bond_dist_bb = None + + if equal_sc: + cgmodel = get_helix_cgmodel_2sc_equal(sigma_bb_init,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) + else: + # Need to set 2 independent particle types for sc1, sc2 + cgmodel = get_helix_cgmodel_2sc_nonequal(sigma_bb_init,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) + + # Get particle type lists and bonded lists: + # (sc1 and sc2 have the same bonded type here) + (particle_type_list, bb_array, sc_array, + bb_bond_list, bs_bond_list, ss_bond_list, + bbb_angle_list, bbs_angle_list, bss_angle_list, + bbbb_torsion_list, bbbs_torsion_list, bbss_torsion_list, + sbbs_torsion_list) = get_helix_particle_bonded_lists_2sc(cgmodel) + + # Set up Simulation object beforehand: + simulation_time_step = 5.0 * unit.femtosecond + friction = 0.0 / unit.picosecond + integrator = LangevinIntegrator( + 0.0 * unit.kelvin, friction, simulation_time_step.in_units_of(unit.picosecond) + ) + simulation = Simulation(cgmodel.topology, cgmodel.system, integrator) + + if bond_dist_bb is not None: + bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) + + if equal_sc: + #-----------------------# + # Equal sidechain sigma # + #-----------------------# + + # Here we also pass sigma_bb to calculate the contact distance for backbone-sc2: + params = (simulation, particle_type_list, r, c, n_particle_bb, bond_dist_bb, sigma_bb, alignment) + + # We can also add the tilt angle defining the cone that sc2 rotates around, + # and make two alternating rotation theta angles + + # Set optimization bounds [sigma_sc, theta1, theta2]: + if n_rotation_angles > 1: + # Allow for n independent rotation of residue sequences + + # Due to symmetry, we can fix the rotation angle bounds to (-pi/2,pi/2). + # This keeps the backbone-sc1 bond on the left side for all residues + + bounds = [] + + if sigma_bb is None: + # Also optimize the backbone parameters: + bounds.append((r/10,10*r)) # sigma_bb + bounds.append((0.01,2*np.pi/3)) # particle spacing in radians + + bounds.append((r/10,10*r)) # sigma_sc + + if alignment == 'center': + for a in range(n_rotation_angles): + bounds.append((-np.pi/2,+np.pi/2)) # Rotation angles + + elif alignment == 'first': + # Due to end-to-end symmetry, we can set one range to span pi + bounds.append((-np.pi/2,+np.pi/2)) + + # The rest of the angles should have full range of rotation: + for a in range(n_rotation_angles-1): + bounds.append((-np.pi,+np.pi)) + + # No constraints if r_ss is fixed by sigma_sc + linear_constraint = () + + # If r_ss is independent variable: + # (Corresponds to [r_ss, sigma_sc, theta1, theta2) + # linear_constraint = LinearConstraint( + # [-1, np.power(2,(1/6)), 0, 0], + # [-np.power(2,(1/6))*sigma_bb.value_in_unit(unit.angstrom)], + # [np.inf] + # ) + + else: + # Uniform rotation across all residues + if sigma_bb is None: + bounds = [(0.01,2*np.pi/3),(r/10,10*r),(r/10,10*r),(-np.pi/2,+np.pi/2)] + # Due to end-to-end symmetry, we can fix the rotation angle bounds to (-pi/2,pi/2) + else: + bounds = [(r/10,10*r),(-np.pi/2,+np.pi/2)] + # No constraints if r_ss is fixed by sigma_sc + linear_constraint = () + + opt_sol = differential_evolution( + compute_helix_openmm_energy_vary_LJ_2sc_rotation_equal, + bounds, + constraints=linear_constraint, + args=params, + polish=True, + popsize=DE_popsize, + ) + + if sigma_bb is None: + t_delta = opt_sol.x[0] + sigma_bb_opt = opt_sol.x[1] * unit.angstrom + + sigma_sc1_opt = opt_sol.x[2] + sigma_sc2_opt = sigma_sc1_opt + + else: + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + sigma_bb_opt = sigma_bb + + sigma_sc1_opt = opt_sol.x[0] + sigma_sc2_opt = sigma_sc1_opt + + theta_opt = [] + + if alignment == 'center' and n_rotation_angles == 2: + + res_per_turn = 2*np.pi/t_delta + n_seq = np.floor(res_per_turn).astype(int) + m_seq = np.ceil(res_per_turn).astype(int) + + if sigma_bb is None: + theta1 = opt_sol.x[3] + theta2 = opt_sol.x[4] + else: + theta1 = opt_sol.x[1] + theta2 = opt_sol.x[2] + i = 0 + while i < n_particle_bb: + for a in range(n_seq): + theta_opt.append(theta1) + i += 1 + for b in range(m_seq): + theta_opt.append(theta2) + i += 1 + + theta_opt = theta_opt[0:n_particle_bb] + + else: + if sigma_bb is None: + for a in range(n_rotation_angles): + theta_opt.append(opt_sol.x[a+3]) + else: + for a in range(n_rotation_angles): + theta_opt.append(opt_sol.x[a+1]) + else: + #---------------------------# + # Non-equal sidechain sigma # + #---------------------------# + # ***TODO: add in the variable backbone option for non-equal sigma_sc + # Here we also pass sigma_bb to calculate the contact distance for backbone-sc2: + params = (simulation, particle_type_list, r, c, n_particle_bb, bond_dist_bb, sigma_bb, alignment) + + # Here we also set r_bs=r_eq_bb_sc1, and the other contact distance will be fixed at r_eq_bb_sc2 + + if n_rotation_angles > 1: + # Set optimization bounds [sigma_sc1, sigma_sc2, theta1, theta2]: + bounds = [] + bounds.append((r/10,10*r)) # sigma_sc1 + bounds.append((r/10,10*r)) # sigma_sc2 + + # For sigma_sc2 >= sigma_sc1: + constraint_var_list = [] + constraint_var_list.append(-1) + constraint_var_list.append(1) + + if alignment == 'center': + for a in range(n_rotation_angles): + bounds.append((-np.pi/2,+np.pi/2)) # Rotation angles + constraint_var_list.append(0) + + elif alignment == 'first': + # Due to end-to-end symmetry, we can set one range to span pi + bounds.append((-np.pi/2,+np.pi/2)) + constraint_var_list.append(0) + + # The rest of the angles should have full range of rotation: + for a in range(n_rotation_angles-1): + bounds.append((-np.pi,+np.pi)) + constraint_var_list.append(0) + + # If we fix r_ss from sigma_sc, there are no constraints: + # linear_constraint = () + + # If we constrain sigma_sc2 >= sigma_sc1, we have the following: + linear_constraint = LinearConstraint( + constraint_var_list, + [0], + [np.inf], + ) + + # If r_ss is independent variable: + # (Corresponds to [r_ss, sigma_sc1, sigma_sc2, theta1, theta2) + # linear_constraint = LinearConstraint( + # [[-1, np.power(2,(-5/6)), np.power(2,(-5/6)), 0, 0], [1, -np.power(2,(-5/6)), np.power(2,(-5/6)), 0, 0], [1, np.power(2,(-5/6)), -np.power(2,(-5/6)), 0, 0]], + # [-np.power(2,(1/6))*sigma_bb.value_in_unit(unit.angstrom), 0, 0], + # [np.inf, np.inf, np.inf] + # ) + else: + # Uniform rotation across all residues: + # Set optimization bounds [sigma_sc1, sigma_sc2, theta1]: + bounds = [(r/10,10*r),(r/10,10*r),(-np.pi/2,+np.pi/2)] + + # If we fix r_ss from sigma_sc, there are no constraints: + # linear_constraint = () + + # If we constrain sigma_sc2 >= sigma_sc1, we have the following: + linear_constraint = LinearConstraint( + [-1, 1, 0], + [0], + [np.inf], + ) + + opt_sol = differential_evolution( + compute_helix_openmm_energy_vary_LJ_2sc_rotation_nonequal, + bounds, + constraints=linear_constraint, + args=params, + polish=True, + popsize=DE_popsize, + ) + + sigma_sc1_opt = opt_sol.x[0] + sigma_sc2_opt = opt_sol.x[1] + + theta_opt = [] + + for a in range(n_rotation_angles): + theta_opt.append(opt_sol.x[a+2]) + + # Compute particle spacing based on bond constraints + if sigma_bb is None: + t_delta_opt = t_delta + else: + t_delta_opt = get_t_from_bond_distance(r,c,bond_dist_bb) + sigma_bb_opt = sigma_bb + + if t_delta_opt < 0: + print(t_delta_opt) + t_delta_opt *= -1 + + t_par = np.zeros(n_particle_bb) + for i in range(n_particle_bb): + t_par[i] = i*t_delta_opt + + # Equilibrium LJ distance (for visual representation) + r_eq_bb = sigma_bb_opt.value_in_unit(unit.angstrom)*np.power(2,(1/6)) + r_eq_sc1 = sigma_sc1_opt*np.power(2,(1/6)) + r_eq_sc2 = sigma_sc2_opt*np.power(2,(1/6)) + + r_eq_bb_sc1 = (sigma_bb_opt.value_in_unit(unit.angstrom)+sigma_sc1_opt)/2*np.power(2,(1/6)) + r_eq_bb_sc2 = (sigma_bb_opt.value_in_unit(unit.angstrom)+sigma_sc2_opt)/2*np.power(2,(1/6)) + r_eq_sc1_sc2 = (sigma_sc1_opt+sigma_sc2_opt)/2*np.power(2,(1/6)) + + # Can set this to r_eq_bb_sc1 for both equal and nonequal sigma_sc cases: + r_bs = r_eq_bb_sc1 + r_ss = r_eq_sc1_sc2 + + # Get particle coorindates: + positions = get_helix_coordinates_2sc_rotation( + r, c, t_par, r_bs, r_ss, r_eq_bb_sc2, theta_opt, alignment + ) + + # Write pdb file + cgmodel.positions = positions + write_pdbfile_without_topology(cgmodel, pdbfile) + + # Also write dcd file (better precision) + dcdfile = pdbfile[:-3]+'dcd' + dcdtraj = md.Trajectory( + xyz=positions, + topology=md.Topology.from_openmm(cgmodel.topology), + ) + md.Trajectory.save_dcd(dcdtraj,dcdfile) + + # Store key geometric parameters + geometry = {} + + geometry['sigma_bb'] = sigma_bb_opt.in_units_of(r_unit) + geometry['sigma_sc1'] = (sigma_sc1_opt*unit.angstrom).in_units_of(r_unit) + geometry['sigma_sc2'] = (sigma_sc2_opt*unit.angstrom).in_units_of(r_unit) + + if alignment == 'center' and n_rotation_angles == 2: + # Use the i-->i+n / i-->i+m sequences + geometry[f'rotation_angle_type1'] = (theta_opt[0]*unit.radian).in_units_of(unit.degrees) + geometry[f'rotation_angle_type2'] = (theta_opt[n_seq]*unit.radian).in_units_of(unit.degrees) + + else: + for a in range(n_rotation_angles): + geometry[f'rotation_angle_type{a+1}'] = (theta_opt[a]*unit.radian).in_units_of(unit.degrees) + + # Add back units: + geometry['helical_radius'] = r * r_unit + geometry['particle_spacing'] = t_delta_opt * unit.radian + geometry['pitch'] = (2*np.pi*c) * r_unit + + # Load dcd file into mdtraj + traj = md.load(dcdfile,top=md.Topology.from_openmm(cgmodel.topology)) + + # Get bb-bb bond distance + geometry['bb_bb_distance'] = (dist_unitless(positions[bb_bond_list[0][0],:],positions[bb_bond_list[0][1],:]) * unit.angstrom).in_units_of(r_unit) + geometry['bb_sc_distance'] = (dist_unitless(positions[bs_bond_list[0][0],:],positions[bs_bond_list[0][1],:]) * unit.angstrom).in_units_of(r_unit) + geometry['sc_sc_distance'] = (dist_unitless(positions[ss_bond_list[0][0],:],positions[ss_bond_list[0][1],:]) * unit.angstrom).in_units_of(r_unit) + + # Get bb-bb-bb angle + angle_indices = np.array([[bbb_angle_list[0][0], bbb_angle_list[0][1], bbb_angle_list[0][2]]]) + geometry['bb_bb_bb_angle'] = (md.compute_angles(traj,angle_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get bb-bb-sc angle + for a in range(n_rotation_angles): + angle_indices = np.array([[bbs_angle_list[a][0], bbs_angle_list[a][1], bbs_angle_list[a][2]]]) + geometry[f'bb_bb_sc_angle_type{a+1}'] = (md.compute_angles(traj,angle_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get bb-sc-sc angle + for a in range(n_rotation_angles): + angle_indices = np.array([[bss_angle_list[a][0], bss_angle_list[a][1], bss_angle_list[a][2]]]) + geometry[f'bb_sc_sc_angle_type{a+1}'] = (md.compute_angles(traj,angle_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get bb-bb-bb-bb torsion + dihedral_indices = np.array([[bbbb_torsion_list[0][0], bbbb_torsion_list[0][1], bbbb_torsion_list[0][2], bbbb_torsion_list[0][3]]]) + geometry['bb_bb_bb_bb_angle'] = (md.compute_dihedrals(traj,dihedral_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get bb-bb-bb-sc torsion + for a in range(n_rotation_angles): + dihedral_indices = np.array([[bbbs_torsion_list[a][0], bbbs_torsion_list[a][1], bbbs_torsion_list[a][2], bbbs_torsion_list[a][3]]]) + geometry[f'bb_bb_bb_sc_angle_type{a+1}'] = (md.compute_dihedrals(traj,dihedral_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get bb-bb-sc-sc torsion + for a in range(n_rotation_angles): + dihedral_indices = np.array([[bbss_torsion_list[a][0], bbss_torsion_list[a][1], bbss_torsion_list[a][2], bbss_torsion_list[a][3]]]) + geometry[f'bb_bb_sc_sc_angle_type{a+1}'] = (md.compute_dihedrals(traj,dihedral_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + + # Get sc-bb-bb-sc torsion (there are n_residues-1 of these) + if n_rotation_angles == n_particle_bb: + for a in range(n_rotation_angles-1): + dihedral_indices = np.array([[sbbs_torsion_list[a][0], sbbs_torsion_list[a][1], sbbs_torsion_list[a][2], sbbs_torsion_list[a][3]]]) + geometry[f'sc_bb_bb_sc_angle_type{a+1}'] = (md.compute_dihedrals(traj,dihedral_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + else: + for a in range(n_rotation_angles): + dihedral_indices = np.array([[sbbs_torsion_list[a][0], sbbs_torsion_list[a][1], sbbs_torsion_list[a][2], sbbs_torsion_list[a][3]]]) + geometry[f'sc_bb_bb_sc_angle_type{a+1}'] = (md.compute_dihedrals(traj,dihedral_indices)*unit.radians).in_units_of(unit.degrees)[0][0] + return opt_sol, geometry + + def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb, sigma_bb, - bond_dist_bb, DE_popsize=200, pdbfile='LJ_helix_3sc_triangle_opt.pdb', - exclusions={}): + bond_dist_bb, DE_popsize=200, pdbfile='LJ_helix_3sc_triangle_opt.pdb', alignment='center', + alternating=True, exclusions={}): """ With specified radius, pitch, backbone sigma, and backbone-backbone bond length, optimize the sidechain sigma and backbone-sidechain bond length in a helix with @@ -775,10 +1315,10 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb :param n_particle_bb: Number of backbone particles to model :type n_particle_bb: int - :param sigma_bb: LJ sigma parameter for backbone beads + :param sigma_bb: LJ sigma parameter for backbone beads. If None, sigma_bb will also be optimized. :type sigma_bb: Quantity - :param bond_dist_bb: bond distance for bb-bb bonds. + :param bond_dist_bb: bond distance for bb-bb bonds. If None, sigma_bb will also be optimized. :type bond_dist_bb: Quantity :param DE_popsize: population size to use in SciPy differential_evolution solver (default=50) @@ -787,6 +1327,12 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb :param pdbfile: Path to pdb file for saving the helical structure (default='LJ_helix_3sc_triangle_opt.pdb') :type pdbfile: str + :param alignment: sidechain alignment scheme - can be 'center' (center of triangle is fixed normal to backbone) or 'first' (first bead is normal to backbone) (default='center') + :type alignment: str + + :param alternating: Option to treat odd and even triangle in-plane rotation independently (default=True) + :type alternating: str + :param exclusions: pass cg_openmm exclusion rules to the cgmodel (by default [0,0,1] is applied to bb-bb, [0,1,1] to bb-sc, sc-sc) :type exclusions: dict @@ -795,6 +1341,11 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb - geometry - Dictionary containing key geometric parameters of the optimized helix """ + # Check input: + if alignment not in ['center','first']: + print(f'Error: invalid alignment input {alignment}') + exit() + r_unit = radius.unit # Use angstrom for writing pdb file: radius = radius.value_in_unit(unit.angstrom) @@ -812,7 +1363,10 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb # Set initial sigma parameters sigma_sc = 1.0 * unit.angstrom - cgmodel = get_helix_cgmodel_triangle(sigma_bb,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) + if sigma_bb is not None: + cgmodel = get_helix_cgmodel_triangle(sigma_bb,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) + else: + cgmodel = get_helix_cgmodel_triangle(sigma_sc,sigma_sc,epsilon_bb,epsilon_sc,n_particle_bb,exclusions) # Get particle type lists and bonded lists: (particle_type_list, bb_array, sc_array, @@ -829,13 +1383,67 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb ) simulation = Simulation(cgmodel.topology, cgmodel.system, integrator) - bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) + if bond_dist_bb is not None: + bond_dist_bb = bond_dist_bb.value_in_unit(unit.angstrom) - params = (simulation, particle_type_list, r, c, n_particle_bb, bond_dist_bb) + params = (simulation, particle_type_list, r, c, n_particle_bb, bond_dist_bb, sigma_bb, alignment, alternating) + + # Set optimization bounds [sigma_sc, theta1, theta2]: - # Set optimization bounds [r_bs, sigma_sc, theta1, theta2]: - bounds = [(r/50,5*r),(r/50,5*r),(0,2*np.pi/3),(0,2*np.pi/3)] + if sigma_bb is not None and bond_dist_bb is not None: + if alternating: + bounds = [(r/50,5*r),(0,2*np.pi/3),(0,2*np.pi/3)] + + # We have a linear constraint that r_ss <= 2*r_bs + if alignment == 'center': + # No constraint needed - only that sigma_sc be positive: + linear_constraint = LinearConstraint( + [1, 0, 0], + [0], + [np.inf] + ) + elif alignment == 'first': + linear_constraint = LinearConstraint( + [(np.power(2,(-1/3))-np.power(2,(1/6))), 0, 0], + [-np.power(2,(-1/3))*sigma_bb.value_in_unit(unit.angstrom)], + [np.inf] + ) + + else: + bounds = [(r/50,5*r),(0,2*np.pi/3)] + + if alignment == 'center': + # No constraint needed - only that sigma_sc be positive: + linear_constraint = LinearConstraint( + [1, 0], + [0], + [np.inf] + ) + elif alignment == 'first': + linear_constraint = LinearConstraint( + [(np.power(2,(-1/3))-np.power(2,(1/6))), 0], + [-np.power(2,(-1/3))*sigma_bb.value_in_unit(unit.angstrom)], + [np.inf] + ) + + elif sigma_bb is None and bond_dist_bb is None: + if alternating: + # [sigma_sc, theta1, theta2, sigma_bb, t] + bounds = [(r/50,5*r),(0,2*np.pi/3),(0,2*np.pi/3),(r/50,15*r),(0.01,2*np.pi/3)] + # TODO: update the constraints for alignment == 'center' + linear_constraint = () + + else: + # [sigma_sc, theta1, sigma_bb, t] + bounds = [(r/50,5*r),(0,2*np.pi/3),(r/50,15*r),(0.01,2*np.pi/3)] + # TODO: update the constraints for alignment == 'center' + linear_constraint = () + + else: + print('Error: sigma_bb and bond_dist_bb must either both be specified, or both be None') + exit() + opt_sol = differential_evolution( compute_helix_openmm_energy_vary_LJ_triangle, bounds, @@ -844,150 +1452,53 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb popsize=DE_popsize, ) - r_bs = opt_sol.x[0] - sigma_sc_opt = opt_sol.x[1] + #r_bs = opt_sol.x[0] + + sigma_sc_opt = opt_sol.x[0] r_ss = sigma_sc_opt*np.power(2,(1/6)) # Angles of rotation in-plane (x axis) for triangle templates 1 and 2 - theta1 = opt_sol.x[2] - theta2 = opt_sol.x[3] + theta1_opt = opt_sol.x[1] + + if alternating: + theta2_opt = opt_sol.x[2] + else: + theta2_opt = theta1_opt - # Compute particle spacing based on bond constraints - t_delta_opt = get_t_from_bond_distance(r,c,bond_dist_bb) + if bond_dist_bb is None: + if alternating: + sigma_bb = opt_sol.x[3] * unit.angstrom + t_delta_opt = opt_sol.x[4] + else: + sigma_bb = opt_sol.x[2] * unit.angstrom + t_delta_opt = opt_sol.x[3] + + else: + # Compute particle spacing based on bond constraints + t_delta_opt = get_t_from_bond_distance(r,c,bond_dist_bb) + if t_delta_opt < 0: print(t_delta_opt) t_delta_opt *= -1 - - t_par = np.zeros(n_particle_bb) + + t = np.zeros(n_particle_bb) for i in range(n_particle_bb): - t_par[i] = i*t_delta_opt + t[i] = i*t_delta_opt # Equilibrium LJ distance (for visual representation) r_eq_bb = sigma_bb.value_in_unit(unit.angstrom)*np.power(2,(1/6)) r_eq_sc = sigma_sc_opt*np.power(2,(1/6)) + r_eq_bb_sc = (sigma_bb.value_in_unit(unit.angstrom)+sigma_sc_opt)/2*np.power(2,(1/6)) - # Get particle positions: - xyz = get_helix_coordinates(r,c,t_par) - - # Place sidechain triangle normal to helix - positions = np.zeros((4*n_particle_bb,3)) - - # From the law of cosines: - if r_bs > r_ss/2: - K = np.sqrt(r_bs**2 - (r_ss**2)/4) # Distance from backbone bead to center of triangle plane - else: - K = np.sqrt(r_ss**2 - (r_bs**2)/4) - - L = np.sqrt((r_ss**2)/3) # Distance from triangle center to lower bead (orientation 1) - M = np.sqrt(L**2 - (r_ss**2)/4) # Distance from triangle center to top of triangle (orientation 1) - - # To get the coordinates of the triangle, we need to rotate the coordinates - # by the angle separating each residue. - - # distance between backbone beads projected onto a circle: - dist_bb_xy = np.sqrt(np.sum(np.power((xyz[0,0:2]-xyz[1,0:2]),2))) - - # helical angle between the two backbone beads projected onto a circle: - theta_arc = np.arccos(1-dist_bb_xy**2/(2*(r**2))) # radians - - # loop over all residues, rotating the reference triangles about the z axis by theta_arc: - - # Orientation 1: - tri_ref_orient1 = np.zeros((3,3)) - - tri_ref_orient1[0,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[0,1] = (1+K/r)*xyz[0,1] - r_ss/2 - tri_ref_orient1[0,2] = xyz[0,2] + M + r_bs = r_eq_bb_sc - tri_ref_orient1[1,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[1,1] = (1+K/r)*xyz[0,1] + r_ss/2 - tri_ref_orient1[1,2] = xyz[0,2] + M - - tri_ref_orient1[2,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[2,1] = (1+K/r)*xyz[0,1] - tri_ref_orient1[2,2] = xyz[0,2] - L - - # Apply the rotation in x: - tri_ref_no_rotate1 = tri_ref_orient1 - - tri_ref_orient1[0,:] = rotate_coordinates_x(tri_ref_orient1[0,:],theta1) - tri_ref_orient1[1,:] = rotate_coordinates_x(tri_ref_orient1[1,:],theta1) - tri_ref_orient1[2,:] = rotate_coordinates_x(tri_ref_orient1[2,:],theta1) - - # Orientation 2: - tri_ref_orient2 = np.zeros((3,3)) - - tri_ref_orient2[0,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[0,1] = (1+K/r)*xyz[0,1] - r_ss/2 - tri_ref_orient2[0,2] = xyz[0,2] - M - - tri_ref_orient2[1,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[1,1] = (1+K/r)*xyz[0,1] + r_ss/2 - tri_ref_orient2[1,2] = xyz[0,2] - M + # Get particle positions: + positions = get_helix_coordinates_3sc_triangle( + r, c, t, r_bs, r_ss, r_eq_bb_sc, theta1_opt, theta2_opt, alignment + ) - tri_ref_orient2[2,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[2,1] = (1+K/r)*xyz[0,1] - tri_ref_orient2[2,2] = xyz[0,2] + L - - # Apply the rotation in x: - tri_ref_no_rotate2 = tri_ref_orient2 - - tri_ref_orient2[0,:] = rotate_coordinates_x(tri_ref_orient2[0,:],theta2) - tri_ref_orient2[1,:] = rotate_coordinates_x(tri_ref_orient2[1,:],theta2) - tri_ref_orient2[2,:] = rotate_coordinates_x(tri_ref_orient2[2,:],theta2) - - z_rise = xyz[1,2] - xyz[0,2] - - j = -1 - for i in range(n_particle_bb): - if i % 2 == 0: - # Orientation 1: - j += 1 - positions[j] = xyz[i] - - j += 1 - triangle_xyz_a = rotate_coordinates_z(tri_ref_orient1[0,:],theta_arc*i) - triangle_xyz_b = rotate_coordinates_z(tri_ref_orient1[1,:],theta_arc*i) - triangle_xyz_c = rotate_coordinates_z(tri_ref_orient1[2,:],theta_arc*i) - - # Use only the x and y from the rotated reference: - positions[j] = triangle_xyz_a - positions[j,2] = triangle_xyz_a[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_b - positions[j,2] = triangle_xyz_b[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_c - positions[j,2] = triangle_xyz_c[2] + z_rise*i - - - else: - # Orientation 2: - j += 1 - positions[j] = xyz[i] - - j += 1 - triangle_xyz_a = rotate_coordinates_z(tri_ref_orient2[0,:],theta_arc*i) - triangle_xyz_b = rotate_coordinates_z(tri_ref_orient2[1,:],theta_arc*i) - triangle_xyz_c = rotate_coordinates_z(tri_ref_orient2[2,:],theta_arc*i) - - # Use only the x and y from the rotated reference: - positions[j] = triangle_xyz_a - positions[j,2] = triangle_xyz_a[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_b - positions[j,2] = triangle_xyz_b[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_c - positions[j,2] = triangle_xyz_c[2] + z_rise*i - - # Write pdb file - cgmodel.positions = positions * unit.angstrom + cgmodel.positions = positions write_pdbfile_without_topology(cgmodel, pdbfile) # Also write dcd file (better precision) @@ -1009,6 +1520,9 @@ def optimize_helix_LJ_parameters_triangle_sidechain(radius, pitch, n_particle_bb geometry['particle_spacing'] = t_delta_opt * unit.radian geometry['pitch'] = (2*np.pi*c) * r_unit + geometry['rotation_angle_type1'] = (theta1_opt*unit.radian).in_units_of(unit.degrees) + geometry['rotation_angle_type2'] = (theta2_opt*unit.radian).in_units_of(unit.degrees) + # Load dcd file into mdtraj traj = md.load(dcdfile,top=md.Topology.from_openmm(cgmodel.topology)) @@ -1099,7 +1613,7 @@ def compute_helix_openmm_energy_vary_LJ(geo, simulation, bb_array, sc_array, for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -1153,25 +1667,58 @@ def compute_helix_openmm_energy_vary_LJ(geo, simulation, bb_array, sc_array, def compute_helix_openmm_energy_vary_LJ_constrained( - geo, simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc): + geo, simulation, bb_array, sc_array, particle_type_list, r, c, n_particle_bb, bond_dist_bb, bond_dist_sc, sigma_bb_in): """ Internal function for computing openmm energy of Lennard-Jones 12-6 helix """ # Backbone sigma parameter - sigma_bb = geo[0] * unit.angstrom + if sigma_bb_in is None: + sigma_bb = geo[0] * unit.angstrom + else: + sigma_bb = sigma_bb_in # Sidechain sigma parameter - sigma_sc = geo[1] * unit.angstrom - - # Particle spacing (radians) - t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + if sigma_bb_in is None: + sigma_sc = geo[1] * unit.angstrom + else: + sigma_sc = geo[0] * unit.angstrom + + if bond_dist_bb is None and bond_dist_sc is not None: + # variable bond_dist_bb, fixed bond_dist_sc + # Particle spacing (radians): + t_delta = geo[2] + r_bb = get_bond_distance_from_t(r,c,t_delta) + r_bs = bond_dist_sc + + elif bond_dist_sc is None and bond_dist_bb is not None: + # variable bond_dist_sc, fixed_bond_dist_bb + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + r_bb = bond_dist_bb + if sigma_bb_in is None: + r_bs = geo[2] + else: + r_bs = geo[1] + + elif bond_dist_sc == 'LJ' and bond_dist_bb is not None: + # set bb-sc bond from sigmas, fixed bond_dist_bb + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + r_bb = bond_dist_bb + r_bs = (sigma_bb.value_in_unit(unit.angstrom)+sigma_sc.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) + + # TODO: add case of 'LJ' bond_dist_sc and variable bond_dist_bb + + else: + # both bond distances fixed + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + r_bb = bond_dist_bb + r_bs = bond_dist_sc t1 = np.zeros(n_particle_bb) for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -1179,8 +1726,6 @@ def compute_helix_openmm_energy_vary_LJ_constrained( # are zero. # Place sidechain particles normal to helix with same bond length as bb_bb - r_bb = bond_dist_bb - r_bs = bond_dist_sc side_xyz = np.zeros((n_particle_bb,3)) side_xyz[:,0] = (1+r_bs/r)*xyz[:,0] @@ -1194,7 +1739,7 @@ def compute_helix_openmm_energy_vary_LJ_constrained( positions[bb_array] = xyz positions[sc_array] = side_xyz - #positions *= unit.angstrom + positions *= unit.angstrom # Update the nonbonded parameters: for force_index, force in enumerate(simulation.system.getForces()): @@ -1236,7 +1781,7 @@ def compute_energy_diff_1sc(distance_vars, simulation, bb_array, sc_array, parti t1[i] = i*t_delta # Get the backbone coordinates: - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -1424,10 +1969,7 @@ def compute_helix_openmm_energy_vary_LJ_2sc_equal(geo, simulation, r_bs = geo[0] # Sidechain-sidechain bond length: - r_ss = geo[1] - - # Backbone sigma parameter: - + r_ss = geo[1] # Sidechain sigma parameter sigma_sc = geo[2] * unit.angstrom @@ -1439,7 +1981,7 @@ def compute_helix_openmm_energy_vary_LJ_2sc_equal(geo, simulation, for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -1522,7 +2064,7 @@ def compute_helix_openmm_energy_vary_LJ_2sc_nonequal(geo, simulation, for i in range(n_particle_bb): t1[i] = i*t_delta - xyz = get_helix_coordinates(r,c,t1) + xyz = get_helix_backbone_coordinates(r,c,t1) # If the bonds, angles, and backbone torsions are at their equilibrium positions, # then we don't need to update any parameters in the simulation object. Just @@ -1582,174 +2124,245 @@ def compute_helix_openmm_energy_vary_LJ_2sc_nonequal(geo, simulation, return U_helix - -def compute_helix_openmm_energy_vary_LJ_triangle( - geo, simulation, particle_type_list, r, c, n_particle_bb, bond_dist_bb): + +def compute_helix_openmm_energy_vary_LJ_2sc_rotation_equal(geo, simulation, + particle_type_list, r, c, n_particle_bb, bond_dist_bb, fixed_sigma_bb, alignment): """ Internal function for computing openmm energy of Lennard-Jones 12-6 helix + (2sc model, equal sigma_sc for sidechain beads) """ - # Backbone-sidechain bond length: - r_bs = geo[0] - - # Sidechain sigma parameter - sigma_sc = geo[1] * unit.angstrom - - # Angle of rotation in x for each triangle reference: - theta1 = geo[2] - theta2 = geo[3] + if fixed_sigma_bb is None: + t_delta = geo[0] + sigma_bb = geo[1] * unit.angstrom + sigma_sc = geo[2] * unit.angstrom + else: + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + sigma_bb = fixed_sigma_bb + sigma_sc = geo[0] * unit.angstrom - # sidechain-sidechain bond length: - r_ss = sigma_sc.value_in_unit(unit.angstrom)*np.power(2,(1/6)) + # Sidechain-sidechain equilibrium distance + r_ss = np.power(2,(1/6))*sigma_sc.value_in_unit(unit.angstrom) - # Particle spacing (radians) - t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + # Equilibrium contact distance for backbone and sidechain bead 2: + r_eq_bb_sc2 = (sigma_bb.value_in_unit(unit.angstrom)+sigma_sc.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) + r_bs = r_eq_bb_sc2 - t1 = np.zeros(n_particle_bb) + t = np.zeros(n_particle_bb) for i in range(n_particle_bb): - t1[i] = i*t_delta - - xyz = get_helix_coordinates(r,c,t1) - - # If the bonds, angles, and backbone torsions are at their equilibrium positions, - # then we don't need to update any parameters in the simulation object. Just - # the nonbonded energies need to be evaluated. In the cgmodel, all force constants - # are zero. - - # Place sidechain triangle normal to helix - positions = np.zeros((4*n_particle_bb,3)) - - # From the law of cosines: - if r_bs > r_ss/2: - K = np.sqrt(r_bs**2 - (r_ss**2)/4) + t[i] = i*t_delta + + res_per_turn = 2*np.pi/t_delta + n_seq = np.floor(res_per_turn).astype(int) + m_seq = np.ceil(res_per_turn).astype(int) + + # Angle of rotation for sidechain 2 bead about backbone-sc1 axis: + theta_list = [] + + # For center alignment hexagonal packing with n < residues/turn < m + # we should apply angles in the following pattern: + # (type1 * n), (type2 * m), (type1 * n), (type2 * m), ... + + if fixed_sigma_bb is None: + if alignment == 'center' and len(geo) == 5: + theta1 = geo[3] + theta2 = geo[4] + i = 0 + while i < n_particle_bb: + for a in range(n_seq): + theta_list.append(theta1) + i += 1 + for b in range(m_seq): + theta_list.append(theta2) + i += 1 + + theta_list = theta_list[0:n_particle_bb] + + else: + # Regularly repeating sequences: + for a in range(1,len(geo)): + theta_list.append(geo[a]) else: - K = np.sqrt(r_ss**2 - (r_bs**2)/4) - - L = np.sqrt((r_ss**2)/3) - M = np.sqrt(L**2 - (r_ss**2)/4) - - # To get the coordinates of the triangle, we need to rotate the coordinates - # by the angle separating each residue. - - # distance between backbone beads projected onto a circle: - dist_bb_xy = np.sqrt(np.sum(np.power((xyz[0,0:2]-xyz[1,0:2]),2))) - - # helical angle between the two backbone beads projected onto a circle: - theta_arc = np.arccos(1-dist_bb_xy**2/(2*(r**2))) # radians - - # loop over all residues, rotating the reference triangles about the z axis by theta_arc: + if alignment == 'center' and len(geo) == 3: + theta1 = geo[1] + theta2 = geo[2] + i = 0 + while i < n_particle_bb: + for a in range(n_seq): + theta_list.append(theta1) + i += 1 + for b in range(m_seq): + theta_list.append(theta2) + i += 1 + + theta_list = theta_list[0:n_particle_bb] + + else: + # Regularly repeating sequences: + for a in range(1,len(geo)): + theta_list.append(geo[a]) + + # Get particle coorindates: + positions = get_helix_coordinates_2sc_rotation( + r, c, t, r_bs, r_ss, r_eq_bb_sc2, theta_list, alignment + ) - # Orientation 1: - tri_ref_orient1 = np.zeros((3,3)) + # Update the nonbonded parameters: + for force_index, force in enumerate(simulation.system.getForces()): + force_name = force.__class__.__name__ + if force_name == 'NonbondedForce': + for particle_index in range(len(particle_type_list)): + (q,sigma_old,eps) = force.getParticleParameters(particle_index) + + # Only need to change the sigma values here: + if particle_type_list[particle_index] == 'bb': + force.setParticleParameters(particle_index,q,sigma_bb,eps) + else: + force.setParticleParameters(particle_index,q,sigma_sc,eps) + force.updateParametersInContext(simulation.context) + + # Update the positions: + simulation.context.setPositions(positions) + potential_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy() - tri_ref_orient1[0,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[0,1] = (1+K/r)*xyz[0,1] - r_ss/2 - tri_ref_orient1[0,2] = xyz[0,2] + M + U_helix = potential_energy.value_in_unit(unit.kilojoules_per_mole) + + return U_helix - tri_ref_orient1[1,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[1,1] = (1+K/r)*xyz[0,1] + r_ss/2 - tri_ref_orient1[1,2] = xyz[0,2] + M - - tri_ref_orient1[2,0] = (1+K/r)*xyz[0,0] - tri_ref_orient1[2,1] = (1+K/r)*xyz[0,1] - tri_ref_orient1[2,2] = xyz[0,2] - L - # Apply the rotation in x: - tri_ref_no_rotate1 = tri_ref_orient1 +def compute_helix_openmm_energy_vary_LJ_2sc_rotation_nonequal(geo, simulation, + particle_type_list, r, c, n_particle_bb, bond_dist_bb, sigma_bb, alignment): + """ + Internal function for computing openmm energy of Lennard-Jones 12-6 helix + (2sc model, equal sigma_sc for sidechain beads) + """ - tri_ref_orient1[0,:] = rotate_coordinates_x(tri_ref_orient1[0,:],theta1) - tri_ref_orient1[1,:] = rotate_coordinates_x(tri_ref_orient1[1,:],theta1) - tri_ref_orient1[2,:] = rotate_coordinates_x(tri_ref_orient1[2,:],theta1) + # Sidechain sigma parameters + sigma_sc1 = geo[0] * unit.angstrom + sigma_sc2 = geo[1] * unit.angstrom - # Orientation 2: - tri_ref_orient2 = np.zeros((3,3)) + # Equilibrium contact distance for backbone and sidechain beads: + r_eq_bb_sc1 = (sigma_bb.value_in_unit(unit.angstrom)+sigma_sc1.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) + r_eq_bb_sc2 = (sigma_bb.value_in_unit(unit.angstrom)+sigma_sc2.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) + r_eq_sc1_sc2 = (sigma_sc1.value_in_unit(unit.angstrom)+sigma_sc2.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) - tri_ref_orient2[0,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[0,1] = (1+K/r)*xyz[0,1] - r_ss/2 - tri_ref_orient2[0,2] = xyz[0,2] - M + r_bs = r_eq_bb_sc1 + r_ss = r_eq_sc1_sc2 - tri_ref_orient2[1,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[1,1] = (1+K/r)*xyz[0,1] + r_ss/2 - tri_ref_orient2[1,2] = xyz[0,2] - M - - tri_ref_orient2[2,0] = (1+K/r)*xyz[0,0] - tri_ref_orient2[2,1] = (1+K/r)*xyz[0,1] - tri_ref_orient2[2,2] = xyz[0,2] + L + # Angle of rotation for sidechain 2 bead about backbone-sc1 axis: + theta_list = [] + for a in range(2,len(geo)): + theta_list.append(geo[a]) - # Apply the rotation in x: - tri_ref_no_rotate2 = tri_ref_orient2 + # Particle spacing (radians) + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) - tri_ref_orient2[0,:] = rotate_coordinates_x(tri_ref_orient2[0,:],theta2) - tri_ref_orient2[1,:] = rotate_coordinates_x(tri_ref_orient2[1,:],theta2) - tri_ref_orient2[2,:] = rotate_coordinates_x(tri_ref_orient2[2,:],theta2) - - z_rise = xyz[1,2] - xyz[0,2] - - j = -1 + t = np.zeros(n_particle_bb) for i in range(n_particle_bb): - if i % 2 == 0: - # Orientation 1: - j += 1 - positions[j] = xyz[i] - - j += 1 - triangle_xyz_a = rotate_coordinates_z(tri_ref_orient1[0,:],theta_arc*i) - triangle_xyz_b = rotate_coordinates_z(tri_ref_orient1[1,:],theta_arc*i) - triangle_xyz_c = rotate_coordinates_z(tri_ref_orient1[2,:],theta_arc*i) - - # Use only the x and y from the rotated reference: - positions[j] = triangle_xyz_a - positions[j,2] = triangle_xyz_a[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_b - positions[j,2] = triangle_xyz_b[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_c - positions[j,2] = triangle_xyz_c[2] + z_rise*i - - - else: - # Orientation 2: - j += 1 - positions[j] = xyz[i] - - j += 1 - triangle_xyz_a = rotate_coordinates_z(tri_ref_orient2[0,:],theta_arc*i) - triangle_xyz_b = rotate_coordinates_z(tri_ref_orient2[1,:],theta_arc*i) - triangle_xyz_c = rotate_coordinates_z(tri_ref_orient2[2,:],theta_arc*i) - - # Use only the x and y from the rotated reference: - positions[j] = triangle_xyz_a - positions[j,2] = triangle_xyz_a[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_b - positions[j,2] = triangle_xyz_b[2] + z_rise*i - - j += 1 - positions[j] = triangle_xyz_c - positions[j,2] = triangle_xyz_c[2] + z_rise*i - - positions *= unit.angstrom + t[i] = i*t_delta + + # Get particle coorindates: + positions = get_helix_coordinates_2sc_rotation( + r, c, t, r_bs, r_ss, r_eq_bb_sc2, theta_list, alignment + ) # Update the nonbonded parameters: for force_index, force in enumerate(simulation.system.getForces()): force_name = force.__class__.__name__ if force_name == 'NonbondedForce': for particle_index in range(len(particle_type_list)): - - # Only need to change the sigma_sc values here: - if particle_type_list[particle_index] == 'sc': - (q,sigma_old,eps) = force.getParticleParameters(particle_index) + (q,sigma_old,eps) = force.getParticleParameters(particle_index) + + # Only need to change the sigma values here: + if particle_type_list[particle_index] == 'bb': + # This is constant and set when creating the cgmodel + pass + elif particle_type_list[particle_index] == 'sc1': + force.setParticleParameters(particle_index,q,sigma_sc1,eps) + elif particle_type_list[particle_index] == 'sc2': + force.setParticleParameters(particle_index,q,sigma_sc2,eps) - force.setParticleParameters(particle_index,q,sigma_sc,eps) - force.updateParametersInContext(simulation.context) + force.updateParametersInContext(simulation.context) # Update the positions: simulation.context.setPositions(positions) potential_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy() + U_helix = potential_energy.value_in_unit(unit.kilojoules_per_mole) + + return U_helix + + +def compute_helix_openmm_energy_vary_LJ_triangle(geo, simulation, + particle_type_list, r, c, n_particle_bb, bond_dist_bb, sigma_bb, alignment, alternating): + """ + Internal function for computing openmm energy of Lennard-Jones 12-6 helix + """ + + # Backbone-sidechain bond length: + # r_bs = geo[0] + + # Sidechain sigma parameter + sigma_sc = geo[0] * unit.angstrom + + # Angle of rotation in x for each triangle reference: + theta1 = geo[1] + + if alternating: + # Independent odd/even residue in-plane triangle rotation (4D optimization) + theta2 = geo[2] + else: + # All residues get the same in-plane triangle rotation (3D optimization) + theta2 = theta1 + + # sidechain-sidechain bond length: + r_ss = sigma_sc.value_in_unit(unit.angstrom)*np.power(2,(1/6)) + + if sigma_bb is None: + if alternating: + sigma_bb_var = geo[3] + t_delta = geo[4] + else: + sigma_bb_var = geo[2] + t_delta = geo[3] + else: + sigma_bb_var = sigma_bb.value_in_unit(unit.angstrom) + # Particle spacing (radians) + t_delta = get_t_from_bond_distance(r,c,bond_dist_bb) + + # sidechain-backbone equilibrium distance: + r_eq_bb_sc = (sigma_bb_var+sigma_sc.value_in_unit(unit.angstrom))/2*np.power(2,(1/6)) + r_bs = r_eq_bb_sc + + t = np.zeros(n_particle_bb) + for i in range(n_particle_bb): + t[i] = i*t_delta + + # Get particle positions: + positions = get_helix_coordinates_3sc_triangle( + r, c, t, r_bs, r_ss, r_eq_bb_sc, theta1, theta2, alignment + ) + + # Update the nonbonded parameters: + try: + for force_index, force in enumerate(simulation.system.getForces()): + force_name = force.__class__.__name__ + if force_name == 'NonbondedForce': + for particle_index in range(len(particle_type_list)): + + # Only need to change the sigma_sc values here: + if particle_type_list[particle_index] == 'sc': + (q,sigma_old,eps) = force.getParticleParameters(particle_index) + + force.setParticleParameters(particle_index,q,sigma_sc,eps) + force.updateParametersInContext(simulation.context) + + # Update the positions: + simulation.context.setPositions(positions) + potential_energy = simulation.context.getState(getEnergy=True).getPotentialEnergy() + U_helix = potential_energy.value_in_unit(unit.kilojoules_per_mole) + + except: + U_helix = 1E6 return U_helix \ No newline at end of file diff --git a/cg_openmm/utilities/helix_utils.py b/cg_openmm/utilities/helix_utils.py index e4ecad80..5768c25b 100644 --- a/cg_openmm/utilities/helix_utils.py +++ b/cg_openmm/utilities/helix_utils.py @@ -51,40 +51,57 @@ def get_extended_positions_1sc(r_bb, r_bs, n_particle_bb, theta_bbb): positions_ext[j+3,1] = h+r_bs j += 4 - return positions_ext - - -def rotate_coordinates_z(xyz, theta): + return positions_ext + +def rotate_coordinates_x(xyz, theta): """ - Internal function for rotating 3d coordinates by theta (angle between adjacent residues) + Internal function for rotating 3d coordinates by theta (angle to rotate the triangle in-plane) Theta is specified in radians. + xyz is formatted as (particle, dimension) """ R_mat = np.array([ - [np.cos(theta), -np.sin(theta), 0], - [np.sin(theta), np.cos(theta), 0], - [0, 0, 1] + [1, 0, 0], + [0, np.cos(theta), -np.sin(theta)], + [0, np.sin(theta), np.cos(theta)], ]) - xyz_rot = np.matmul(R_mat, xyz) + xyz_rot = np.matmul(R_mat, xyz) return xyz_rot -def rotate_coordinates_x(xyz, theta): +def rotate_coordinates_y(xyz, theta): """ Internal function for rotating 3d coordinates by theta (angle to rotate the triangle in-plane) Theta is specified in radians. """ R_mat = np.array([ - [1, 0, 0], - [0, np.cos(theta), -np.sin(theta)], - [0, np.sin(theta), np.cos(theta)], + [np.cos(theta), 0, np.sin(theta)], + [0, 1, 0], + [-np.sin(theta), 0, np.cos(theta)], ]) xyz_rot = np.matmul(R_mat, xyz) + return xyz_rot + + +def rotate_coordinates_z(xyz, theta): + """ + Internal function for rotating 3d coordinates by theta (angle between adjacent residues) + Theta is specified in radians. + """ + + R_mat = np.array([ + [np.cos(theta), -np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 0], + [0, 0, 1] + ]) + + xyz_rot = np.matmul(R_mat, xyz) + return xyz_rot @@ -192,7 +209,7 @@ def get_helix_cgmodel(sigma_bb, sigma_sc, epsilon_bb, epsilon_sc, n_particle_bb, for i in range(n_particle_bb): t0[i] = i*np.pi/5 - xyz_bb = get_helix_coordinates(1,1,t0) + xyz_bb = get_helix_backbone_coordinates(1,1,t0) xyz_sc = np.zeros((n_particle_bb,3)) @@ -336,7 +353,7 @@ def get_helix_cgmodel_2sc_equal(sigma_bb, sigma_sc, epsilon_bb, epsilon_sc, n_pa for i in range(n_particle_bb): t0[i] = i*np.pi/5 - xyz_bb = get_helix_coordinates(1,1,t0) + xyz_bb = get_helix_backbone_coordinates(1,1,t0) j = -1 for i in range(n_particle_bb): @@ -495,7 +512,7 @@ def get_helix_cgmodel_2sc_nonequal(sigma_bb, sigma_sc, epsilon_bb, epsilon_sc, n for i in range(n_particle_bb): t0[i] = i*np.pi/5 - xyz_bb = get_helix_coordinates(1,1,t0) + xyz_bb = get_helix_backbone_coordinates(1,1,t0) j = -1 for i in range(n_particle_bb): @@ -645,7 +662,7 @@ def get_helix_cgmodel_triangle(sigma_bb, sigma_sc, epsilon_bb, epsilon_sc, n_par for i in range(n_particle_bb): t0[i] = i*np.pi/5 - xyz_bb = get_helix_coordinates(1,1,t0) + xyz_bb = get_helix_backbone_coordinates(1,1,t0) j = -1 for i in range(n_particle_bb): @@ -995,9 +1012,9 @@ def get_helix_particle_bonded_lists_triangle(cgmodel): bsss_torsion_list, sbbs_torsion_list) -def get_helix_coordinates(r, c, t): +def get_helix_backbone_coordinates(r, c, t): """ - Internal functon for getting the coordinates of particles along a helix, + Internal functon for getting the backbone coordinates of particles along a helix, with positions t. """ @@ -1010,6 +1027,391 @@ def get_helix_coordinates(r, c, t): return xyz +def get_helix_coordinates_2sc_rotation(r, c, t, r_bs, r_ss, r_eq_bb_sc2, theta, alignment): + """ + Internal functon for getting the full coordinates of particles along a helix + (2sc model with bent bb-sc1-sc2 angle) + + :param r: helical radius, in Angstrom + :type r: float + + :param c: helical rise parameter, in Angstrom + :type c: float + + :param t: backbone particle spacing, in radians + :type t: float + + :param r_bs: backbone-sidechain bond length, in Angstrom + :type r_bs: float + + :param r_ss: sidechain-sidechain bond length, in Angstrom + :type r_ss: float + + :param r_eq_bb_sc2: equilibrium contact distance for backbone and 2nd sidechain bead, in Angstrom + :type r_eq_bb_sc2: float + + :param theta: angle(s) of rotation for sc2 about bb-sc1 axis, in radians + :type theta: list(float) + + :param alignment: sidechain alignment scheme - can be 'center' (center of sidechain group is fixed normal to backbone) or 'first' (first bead is normal to backbone) + :type alignment: str + + :returns: + - positions - np.array( n_particles_bb*3 x 3 ) + + """ + + n_particle_bb = len(t) + + # Get backbone coordinates + xyz_backbone = get_helix_backbone_coordinates(r,c,t) + + ref_orient = {} + # Type 1 rotation angle: + ref_orient[1] = np.zeros((3,3)) + + # Backbone positions in first residue: + ref_orient[1][0,:] = xyz_backbone[0,:] + + if alignment == 'first': + # Sidechain bead 1 positions in first residue, normal to backbone: + ref_orient[1][1,0] = (1+r_bs/r)*xyz_backbone[0,0] + ref_orient[1][1,1] = (1+r_bs/r)*xyz_backbone[0,1] + ref_orient[1][1,2] = xyz_backbone[0,2] + + # Place second sidechain particle at optimal local rotation: + # This theta_sbs tilt angle could also be set as an optimization variable for more degrees of freedom + if (r_eq_bb_sc2**2 + r_bs**2 - r_ss**2)/(2*r_eq_bb_sc2*r_bs) > 1: + print(f'arccos(x), x > 1') + print(f'r_eq_bb_sc2: {r_eq_bb_sc2}') + print(f'r_bs: {r_bs}') + print(f'r_ss: {r_ss}') + print(f'num: {(r_eq_bb_sc2**2 + r_bs**2 - r_ss**2)}') + print(f'den: {(2*r_eq_bb_sc2*r_bs)}') + + elif (r_eq_bb_sc2**2 + r_bs**2 - r_ss**2)/(2*r_eq_bb_sc2*r_bs) < -1: + print(f'arccos(x), x < -1') + print(f'r_eq_bb_sc2: {r_eq_bb_sc2}') + print(f'r_bs: {r_bs}') + print(f'r_ss: {r_ss}') + print(f'num: {(r_eq_bb_sc2**2 + r_bs**2 - r_ss**2)}') + print(f'den: {(2*r_eq_bb_sc2*r_bs)}') + + theta_sbs = np.arccos((r_eq_bb_sc2**2 + r_bs**2 - r_ss**2)/(2*r_eq_bb_sc2*r_bs)) + + # Shift in x from the first backbone particle: + # This is constant for a tilt angle theta_sbs + # xshift should always be positive + xshift = np.abs(r_eq_bb_sc2*np.cos(theta_sbs)) + + # Shift in y from the first backbone particle: + # This needs to be rotated about the y axis by theta: + yshift = r_eq_bb_sc2*np.sin(theta_sbs) + + ref_orient[1][2,0] = xyz_backbone[0,0] + xshift + ref_orient[1][2,1] = xyz_backbone[0,1] + yshift + ref_orient[1][2,2] = xyz_backbone[0,2] + + elif alignment == 'center': + # Sidechain group center is normal to backbone: + + # xshift should always be positive and real + xshift = np.sqrt(r_bs**2 - (r_ss**2)/4) + yshift = r_ss/2 + + # sc1 particle: + ref_orient[1][1,0] = xyz_backbone[0,0] + xshift + ref_orient[1][1,1] = xyz_backbone[0,1] - yshift + ref_orient[1][1,2] = xyz_backbone[0,2] + + # sc2 particle: + ref_orient[1][2,0] = xyz_backbone[0,0] + xshift + ref_orient[1][2,1] = xyz_backbone[0,1] + yshift + ref_orient[1][2,2] = xyz_backbone[0,2] + + # Here we can have any number of rotation angles: + n_rotation_angles = len(theta) + for a in range(1,n_rotation_angles): + ref_orient[a+1] = ref_orient[1] + + # Apply the rotation about x axis - this only changes the positions of sidechain bead 2: + # The matrix rows here should be x, y, z and columns the particle indices + + ref_orient_rotx = {} + ref_orient_rotx[1] = np.transpose(rotate_coordinates_x(np.transpose(ref_orient[1]),theta[0])) + + for a in range(1,n_rotation_angles): + ref_orient_rotx[a+1] = np.transpose(rotate_coordinates_x(np.transpose(ref_orient[1]),theta[a])) + + # Now, rotate the template residue about the z axis to construct the full helix: + positions = np.zeros((3*n_particle_bb,3)) + + # distance between backbone beads projected onto a circle: + dist_bb_xy = np.sqrt(np.sum(np.power((xyz_backbone[0,0:2]-xyz_backbone[1,0:2]),2))) + + # helical angle between the two backbone beads projected onto a circle: + theta_arc = np.arccos(1-dist_bb_xy**2/(2*(r**2))) # radians + + # Vertical rise from one residue to the next: + z_rise = xyz_backbone[1,2] - xyz_backbone[0,2] + + # Now, assign each residue to a rotation angle index + rotation_ids = [] + # theta_set = np.unique(theta) + + # if len(theta_set) < len(theta): + # # Special angle sequence for i-->i+n, i-->i+m + # for i in range(n_particle_bb): + # # TODO: generalize to n thetas instead of 2 + # if theta[i] == theta_set[0]: + # rotation_ids.append(0) + # elif theta[i] == theta_set[1]: + # rotation_ids.append(1) + + i = 0 + while i < n_particle_bb: + for a in range(n_rotation_angles): + rotation_ids.append(a) + i += 1 + + rotation_ids = rotation_ids[0:n_particle_bb] + + # Now, do the rotation about z axis to construct the helix: + j = -1 + for i in range(n_particle_bb): + ref_rotx_curr = eval(f'ref_orient_rotx[rotation_ids[i]+1]') + + ref_orient_rotx_rotz_curr = np.transpose( + rotate_coordinates_z(np.transpose(ref_rotx_curr),theta_arc*i)) + + # The z coordinate must be also shifted by the rise/residue + + # Backbone + j += 1 + positions[j,:] = ref_orient_rotx_rotz_curr[0,:] + positions[j,2] = ref_orient_rotx_rotz_curr[0,2] + z_rise*i + + # Sidechain 1 + j += 1 + positions[j,:] = ref_orient_rotx_rotz_curr[1,:] + positions[j,2] = ref_orient_rotx_rotz_curr[1,2] + z_rise*i + + # Sidechain 2 + j += 1 + positions[j,:] = ref_orient_rotx_rotz_curr[2,:] + positions[j,2] = ref_orient_rotx_rotz_curr[2,2] + z_rise*i + + positions *= unit.angstrom + + return positions + + +def get_helix_coordinates_3sc_triangle(r, c, t, r_bs, r_ss, r_eq_bb_sc, theta1, theta2, alignment): + """ + Internal functon for getting the full coordinates of particles along a helix + (2sc model with bent bb-sc1-sc2 angle) + + :param r: helical radius, in Angstrom + :type r: float + + :param c: helical rise parameter, in Angstrom + :type c: float + + :param t: backbone particle spacing, in radians + :type t: float + + :param r_bs: backbone-sidechain bond length, in Angstrom + :type r_bs: float + + :param r_ss: sidechain-sidechain bond length, in Angstrom + :type r_ss: float + + :param r_eq_bb_sc: equilibrium contact distance for backbone and non-bonded sidechain bead, in Angstrom + :type r_eq_bb_sc: float + + :param theta1: angle of rotation for sc2 about bb-sc1 axis (even residues), in radians + :type theta1: float + + :param theta2: angle of rotation for sc2 about bb-sc1 axis (odd residues), in radians + :type theta2: float + + :param alignment: sidechain alignment scheme - can be 'center' (center of sidechain group is fixed normal to backbone) or 'first' (first bead is normal to backbone) + :type alignment: str + + :returns: + - positions - np.array( n_particles_bb*3 x 3 ) + + """ + n_particle_bb = len(t) + + # Get backbone coordinates + xyz_backbone = get_helix_backbone_coordinates(r,c,t) + + ref_orient1 = np.zeros((4,3)) + + # Backbone positions in first residue: + ref_orient1[0,:] = xyz_backbone[0,:] + + if alignment == 'first': + # Sidechain bead 1 positions in first residue, normal to backbone: + ref_orient1[1,0] = xyz_backbone[0,0] + ref_orient1[1,1] = xyz_backbone[0,1] + r_bs + ref_orient1[1,2] = xyz_backbone[0,2] + + # Due to shift in z for sidechains 2,3, need the distance between backbone + # and center of line connecting sc2-sc3: + + d_bb_sc23 = np.sqrt(r_eq_bb_sc**2-(r_ss**2)/4) + + # r_ss distance projected onto the xy plane: + r_ss_xy = r_ss*np.sqrt(3)/2 + + if (d_bb_sc23**2 + r_bs**2 - r_ss_xy**2)/(2*d_bb_sc23*r_bs) > 1: + print(f'arccos(x), x > 1') + print(f'd_bb_sc23: {d_bb_sc23}') + print(f'r_bs: {r_bs}') + print(f'r_ss: {r_ss}') + print(f'r_ss_xy: {r_ss_xy}') + print(f'num: {(d_bb_sc23**2 + r_bs**2 - r_ss_xy**2)}') + print(f'den: {(2*d_bb_sc23*r_bs)}') + + elif (d_bb_sc23**2 + r_bs**2 - r_ss_xy**2)/(2*d_bb_sc23*r_bs) < -1: + print(f'arccos(x), x < -1') + print(f'd_bb_sc23: {d_bb_sc23}') + print(f'r_bs: {r_bs}') + print(f'r_ss: {r_ss}') + print(f'r_ss_xy: {r_ss_xy}') + print(f'num: {(d_bb_sc23**2 + r_bs**2 - r_ss_xy**2)}') + print(f'den: {(2*d_bb_sc23*r_bs)}') + + # Angle in xy plane for sidechain-backbone-sidechain + # theta_sbs = np.arccos((d_bb_sc23**2 + r_bs**2 - r_ss_xy**2)/(2*d_bb_sc23*r_bs)) + + # Shift in x from the first backbone particle: + # This is constant for a tilt angle theta_sbs + # xshift should always be positive + xshift = np.sqrt(3/4*r_ss**2-1/4*r_ss**4/r_bs**2) + + # Shift in y from the first backbone particle: + # This needs to be rotated about the y axis by theta: + yshift = r_bs - 1/2*r_ss**2/r_bs + + # For equilateral triangle model, also shift in z by r_ss/2: + zshift = r_ss/2 + + # sc2 particle: + ref_orient1[2,0] = xyz_backbone[0,0] + xshift + ref_orient1[2,1] = xyz_backbone[0,1] + yshift + ref_orient1[2,2] = xyz_backbone[0,2] + zshift + + # sc3 particle: + ref_orient1[3,0] = xyz_backbone[0,0] + xshift + ref_orient1[3,1] = xyz_backbone[0,1] + yshift + ref_orient1[3,2] = xyz_backbone[0,2] - zshift + + + elif alignment == 'center': + # Sidechain group center is normal to backbone: + + K = np.sqrt(r_bs**2 - (r_ss**2)/4) # Distance from backbone bead to midpoint of sc1-sc2 bond + L = np.sqrt((r_ss**2)/3) # Distance from triangle center to lower bead (sc3) + M = np.sqrt((r_ss**2)/12) # Distance from triangle center to top of triangle (sc1, sc2 beads) + + xshift = np.sqrt(K**2 - M**2) + + # sc1 particle: + ref_orient1[1,0] = xyz_backbone[0,0] + xshift + ref_orient1[1,1] = xyz_backbone[0,1] - r_ss/2 + ref_orient1[1,2] = xyz_backbone[0,2] + M + + # sc2 particle: + ref_orient1[2,0] = xyz_backbone[0,0] + xshift + ref_orient1[2,1] = xyz_backbone[0,1] + r_ss/2 + ref_orient1[2,2] = xyz_backbone[0,2] + M + + # sc3 particle: + ref_orient1[3,0] = xyz_backbone[0,0] + xshift + ref_orient1[3,1] = xyz_backbone[0,1] + ref_orient1[3,2] = xyz_backbone[0,2] - L + + ref_orient2 = ref_orient1 + + # Apply the rotation about x axis - this only changes the positions of sidechain bead 2: + # The matrix rows here should be x, y, z and columns the particle indices + + ref_orient1_rotx = np.transpose(rotate_coordinates_x(np.transpose(ref_orient1),theta1)) + ref_orient2_rotx = np.transpose(rotate_coordinates_x(np.transpose(ref_orient2),theta2)) + + # Now, rotate the template residue about the z axis to construct the full helix: + positions = np.zeros((4*n_particle_bb,3)) + + # distance between backbone beads projected onto a circle: + dist_bb_xy = np.sqrt(np.sum(np.power((xyz_backbone[0,0:2]-xyz_backbone[1,0:2]),2))) + + # helical angle between the two backbone beads projected onto a circle: + theta_arc = np.arccos(1-dist_bb_xy**2/(2*(r**2))) # radians + + # Vertical rise from one residue to the next: + z_rise = xyz_backbone[1,2] - xyz_backbone[0,2] + + j = -1 + for i in range(n_particle_bb): + if i % 2 == 0: + # Orientation 1 + ref_orient1_rotx_rotz = np.transpose(rotate_coordinates_z(np.transpose(ref_orient1_rotx),theta_arc*i)) + + # The z coordinate must be also shifted by the rise/residue + + # Backbone + j += 1 + positions[j,:] = ref_orient1_rotx_rotz[0,:] + positions[j,2] = ref_orient1_rotx_rotz[0,2] + z_rise*i + + # Sidechain 1 + j += 1 + positions[j,:] = ref_orient1_rotx_rotz[1,:] + positions[j,2] = ref_orient1_rotx_rotz[1,2] + z_rise*i + + # Sidechain 2 + j += 1 + positions[j,:] = ref_orient1_rotx_rotz[2,:] + positions[j,2] = ref_orient1_rotx_rotz[2,2] + z_rise*i + + # Sidechain 3 + j += 1 + positions[j,:] = ref_orient1_rotx_rotz[3,:] + positions[j,2] = ref_orient1_rotx_rotz[3,2] + z_rise*i + + else: + # Orientation 2 + ref_orient2_rotx_rotz = np.transpose(rotate_coordinates_z(np.transpose(ref_orient2_rotx),theta_arc*i)) + + # The z coordinate must be also shifted by the rise/residue + + # Backbone + j += 1 + positions[j,:] = ref_orient2_rotx_rotz[0,:] + positions[j,2] = ref_orient2_rotx_rotz[0,2] + z_rise*i + + # Sidechain 1 + j += 1 + positions[j,:] = ref_orient2_rotx_rotz[1,:] + positions[j,2] = ref_orient2_rotx_rotz[1,2] + z_rise*i + + # Sidechain 2 + j += 1 + positions[j,:] = ref_orient2_rotx_rotz[2,:] + positions[j,2] = ref_orient2_rotx_rotz[2,2] + z_rise*i + + # Sidechain 3 + j += 1 + positions[j,:] = ref_orient2_rotx_rotz[3,:] + positions[j,2] = ref_orient2_rotx_rotz[3,2] + z_rise*i + + positions *= unit.angstrom + + return positions + def get_t_from_bond_distance(r, c, bond_dist_bb): """ Internal function for calculating the change in t (from helix parameter equation) @@ -1063,10 +1465,10 @@ def plot_LJ_helix(r, c, t_par, r_eq_bb, r_eq_sc=None, plotfile='LJ_helix.pdf'): # TODO: add sidechain particles to the 3d plot # Helix coordinates: - xyz_helix = get_helix_coordinates(r,c,np.linspace(0,1.5*max(t_par),num=1001)) + xyz_helix = get_helix_backbone_coordinates(r,c,np.linspace(0,1.5*max(t_par),num=1001)) # Particle coordinates: - xyz_par = get_helix_coordinates(r,c,t_par) + xyz_par = get_helix_backbone_coordinates(r,c,t_par) fig = plt.figure()