def run_simulation(force_constant):
    assert(os.path.exists(folder_to_store_output_files))
    input_pdb_file_of_molecule = args.starting_pdb_file
    force_field_file = 'amber99sb.xml'
    water_field_file = 'tip3p.xml'
    pdb_reporter_file = '%s/output_fc_%f_pc_%s.pdb' %(folder_to_store_output_files, force_constant, str(potential_center).replace(' ',''))

    if not args.out_traj is None:
        pdb_reporter_file = args.out_traj

    state_data_reporter_file = pdb_reporter_file.replace('output_fc', 'report_fc').replace('.pdb', '.txt')

    # check if the file exist
    for item_filename in [pdb_reporter_file, state_data_reporter_file]:
        Helper_func.backup_rename_file_if_exists(item_filename)

    index_of_backbone_atoms = CONFIG_57[0]
    flag_random_seed = 0 # whether we need to fix this random seed

    simulation_temperature = args.temperature
    time_step = CONFIG_22   # simulation time step, in ps

    pdb = PDBFile(input_pdb_file_of_molecule)
    modeller = Modeller(pdb.topology, pdb.getPositions(frame=args.starting_frame))
    solvent_opt = 'no_water'
    if solvent_opt == 'explicit':
        forcefield = ForceField(force_field_file, water_field_file)
        modeller.addSolvent(forcefield, model=water_field_file.split('.xml')[0], boxSize=Vec3(3, 3, 3) * nanometers,
                            ionicStrength=0 * molar)
        system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=1.0 * nanometers,
                                         constraints=AllBonds, ewaldErrorTolerance=0.0005)
    else:
        forcefield = ForceField(force_field_file)
        system = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff, constraints=AllBonds)

    if args.bias_method == "US":
        if float(force_constant) != 0:
            force = ANN_Force()
            force.set_layer_types(layer_types)
            force.set_data_type_in_input_layer(args.data_type_in_input_layer)
            force.set_list_of_index_of_atoms_forming_dihedrals_from_index_of_backbone_atoms(index_of_backbone_atoms)
            force.set_index_of_backbone_atoms(index_of_backbone_atoms)
            if args.data_type_in_input_layer == 2:
                force.set_list_of_pair_index_for_distances(CONFIG_80)

            force.set_num_of_nodes(num_of_nodes)
            force.set_potential_center(potential_center)
            force.set_force_constant(float(force_constant))
            unit_scaling = 1.0  # TODO: check unit scaling
            force.set_scaling_factor(float(scaling_factor) / unit_scaling)   # since default unit is nm in OpenMM

            # TODO: need to fix following for multi-hidden layer cases
            temp_coeffs, temp_bias = np.load(autoencoder_info_file)
            for item_layer_index in [0, 1]:
                assert (len(temp_coeffs[item_layer_index]) ==
                        num_of_nodes[item_layer_index] * num_of_nodes[item_layer_index + 1]), (len(temp_coeffs[item_layer_index]),
                                (num_of_nodes[item_layer_index], num_of_nodes[item_layer_index + 1]))
                assert (len(temp_bias[item_layer_index]) == num_of_nodes[item_layer_index + 1]), (len(temp_bias[item_layer_index]), num_of_nodes[item_layer_index + 1])

            # need tolist() since C++ only accepts Python list
            force.set_coeffients_of_connections([item_w.tolist() for item_w in temp_coeffs])
            force.set_values_of_biased_nodes([item_w.tolist() for item_w in temp_bias])

            system.addForce(force)
    elif args.bias_method == "US_on_phipsi":
        from openmmplumed import PlumedForce
        kappa_string = ','.join([str(force_constant) for _ in potential_center])
        plumed_force_string = """