Are restraints removed?

[pipitoludovico]

OpenMM version: 8.0.0 Python V: 3.10

Hi! A quick check: given this equilibration

`def Equilibrate(simulation_eq, parameters) -> None: protocolSteps = int((parameters['protocolTime'] * 10 * 6) / parameters['eqIntegrator']) print('\nEquilibrating with restraints and temperature for ', protocolSteps, " steps.") equilibrationSteps = int((parameters['equilibrationTime'] 10 ** 6) / parameters['eqIntegrator'])

easing = int(protocolSteps * 0.2)
prot = int(protocolSteps - easing)
simulation_eq.step(prot)
print(f'\n{easing} steps easing "k" restraint...:')
count = 0
while count < easing:
    for i_T in range(50, -1, -2):
        r_i = (round(i_T * 0.1, 3))
        r_i = max(r_i, 0)
        print('restraint coefficient = ', r_i)
        simulation_eq.context.setParameter('k', r_i * kilocalories_per_mole / nanometer ** 2)
        simulation_eq.step(1000)
        count += 1000

simulation_eq.context.setParameter('k', 0)
simulation_eq.step(equilibrationSteps)
final_state = simulation_eq.context.getState(getPositions=True, getVelocities=True, enforcePeriodicBox=True)
with open('equilibration_checkpnt.xml', 'w') as output:
    output.write(XmlSerializer.serialize(final_state))
simulation_eq.reporters.clear()`

I want to run the production with no restraints at all. The xml file should have no global parameters saved and therefore when called the new simulation should be restraint free. Is this correct?

def RunProduction(GPU, replica, system, top, parameters):
    productionSteps = int((parameters['productionTime'] * 10 ** 6) / parameters['prIntegrator'])
    integrator, platform, properties = GetMDsettings()
    properties = {'DeviceIndex': f'{str(GPU)}', 'Precision': 'mixed'}
    makedirs(f"run_{replica}", exist_ok=True)
    simulation_production = Simulation(top.topology, system, integrator, platform, properties)
    if os.path.exists("equilibration_checkpnt.xml"):
        simulation_production.loadState("equilibration_checkpnt.xml")
    else:
        raise FileNotFoundError('No XML file found. Make sure the XmlSerializer wrote the state file.')

    simulation_production.reporters.append(CheckpointReporter(f'run_{replica}/{replica}_checkpnt.chk', parameters['saveFreq']))
    simulation_production.reporters.append(DCDReporter(f"run_{replica}/Traj_{replica}.dcd", parameters['saveFreq'], enforcePeriodicBox=True))
    simulation_production.reporters.append(StateDataReporter(f'run_{replica}/Traj_{replica}.std', 1000, step=True, totalSteps=productionSteps, remainingTime=True, potentialEnergy=True, temperature=True, speed=True))

    simulation_production.step(productionSteps)

    final_state = simulation_production.context.getState(getPositions=True, getVelocities=True, enforcePeriodicBox=True)
    with open(f'run_{replica}/Traj_{replica}.xml', 'w') as output:
        output.write(XmlSerializer.serialize(final_state))
    simulation_production.reporters.clear()

Ideally it should have the global parameter 'k' inherited to 0, but it seems that the RMSD of my system is very low and I'm afraid it has not removed the restraints at all.

Thanks

[peastman]

Since you didn't include your complete code, I can only guess at it. I gather you have some sort of a restraint force that defines a global parameter k? If you want the saved state to include the values of parameters, you should specify getParameters=True in the call to getState(). Since you didn't, the state doesn't include parameters and the call to loadState() doesn't change them. The parameter k still has whatever the default value was that you specified when you created the restraint force.

If you're ever in doubt about what the value of a parameter is, you can just print out simulation.context.getParameter('k').

[pipitoludovico]

Everything sorted. Thank you!