Dear Alain,

Thanks for your message and for suggesting the fix for the PULCHRA installation – I've updated the notebook accordingly.
We've also recently made available a new version of the Colab notebook, which implements simulations of IDPs and multi-domain proteins using the latest version of the model.

Regarding the restraint between residues 1 and 66, you could apply it as a harmonic potential, adding the lines below in the the simulate function within the MD Toolbox cell of the Colab notebook.

  hb_cc = openmm.openmm.HarmonicBondForce()
  res_i = 0
  res_j = 65
  r_eq = 0.4*unit.nanometer
  force_constant = 1e4*unit.kilojoules_per_mole/(unit.nanometer**2) # ballpark value!
  hb_cc.addBond(res_i, res_j, r_eq, force_constant)
  yu.addExclusion(res_i, res_j)
  ah.addExclusion(res_i, res_j)
  hb_cc.setUsesPeriodicBoundaryConditions(True)

I hope this helps!

Best regards,
Giulio

Hello Giulio

Many thanks for your update and the new collab. I ll test asap on our 'protein' of interest where the folding process is quite well known and ordered.
That will help us a lot

add : test on first colab with constraints, was fine.
I ll try soon on the new colab, using a pdb of the folded domain and a constraint.

Best regards

Alain

.....
hb = openmm.openmm.HarmonicBondForce()
energy_expression = 'select(step(r-2^(1/6)s),4epsl((s/r)^12-(s/r)^6-shift),4eps((s/r)^12-(s/r)^6-lshift)+eps(1-l))'
ah = openmm.openmm.CustomNonbondedForce(energy_expression+'; s=0.5*(s1+s2); l=0.5*(l1+l2); shift=(0.5*(s1+s2)/2)^12-(0.5*(s1+s2)/2)^6')
yu = openmm.openmm.CustomNonbondedForce('q*(exp(-kappar)/r - exp(-kappa4)/4); q=q1*q2')
yu.addGlobalParameter('kappa',yukawa_kappa/unit.nanometer)
yu.addPerParticleParameter('q')

hb_cc = openmm.openmm.HarmonicBondForce()
res_i = 5
res_j = 61
r_eq = 0.4*unit.nanometer
force_constant = 1e4*unit.kilojoules_per_mole/(unit.nanometer**2) # ballpark value!
hb_cc.addBond(res_i, res_j, r_eq, force_constant)
yu.addExclusion(res_i, res_j)
ah.addExclusion(res_i, res_j)
hb_cc.setUsesPeriodicBoundaryConditions(True)

ah.addGlobalParameter('eps',lj_eps*unit.kilojoules_per_mole)
ah.addPerParticleParameter('s')
ah.addPerParticleParameter('l')

for a,e in zip(seq,yukawa_eps):