[WIP] Extended dynamics for harmonic restraints

[fhh2626]

As mentioned in #649 , this PR makes it possible to use a harmonic restraint as a collective variable (CV). Enhanced sampling algorithms can then be used to calculate the free-energy change associated with creating or removing the harmonic restraint(s).

Examples:

colvar {
  name d
  distance {
    group1 {
      atomnumbers   4759
    }
    group2 {
      dummyAtom (-2.543, 1.510, 0.364)    
    }
  }
}
harmonic {
  name       harm
  colvars       d
  centers       0.0
  forceConstant 10.0
}

colvar {
  name k_cv

  extendedLagrangian on
  extendedMass 1500000
  extendedLangevinDamping 5000

  lowerBoundary      0.0
  upperBoundary      1.0
  width 0.01
  
  reflectinglowerboundary  on
  reflectingupperboundary  on

  harmonicForceConstant {
    harmonicName harm
    kExponent    4.0
  }
}

 abf {
   colvars k_cv
   fullSamples 1000
   outputFreq  500000
   historyFreq 500000
}
metadynamics {
    colvars             k_cv
    hillWeight          0.05
    hillWidth           3.0
    wellTempered        on
    biasTemperature     4000
    keepFreeEnergyFiles on
    outputFreq   500000
}

Next steps:

• Extensive tests

[jhenin]

Thanks @fhh2626 ! I don't have time for an in-depth review right now, but as a reminder to myself I'll write here some points I'd like to watch for:

• maximum code reuse with the alchemical coordinate
• compatibility with the recently overhauled moving restraint code
• extensibility to other parameters e.g. harmonic wall position

[jhenin]

I think it should be possible to avoid these intrusive changes to the colvar class, if the linkage is driven by the bias instead. That is, at init time, the bias could be given a colvar to link to, look it up by name (that is already implemented) and set the mutual pointers.

[fhh2626]

Thanks for your comment. I've removed the link_biases function and done the reverse lookup, but I am not sure how to bypass get_cvc_ptr.

[jhenin]

Unless I'm misunderstanding this, the force constant update should be done once, not once per colvar, so it should not be inside this loop.

More generally, this reimplements colvarbias_restraint::update() with a lot of duplication. Instead, this should be implemented as a variant of the moving restraint code: this is another way of updating the force constant. It would be very clear if the update was done in a similar way, calling the linked cv instead of computing the new force constant internally.

[fhh2626]

I have reused the current implementation of update().

[fhh2626]

Please have a look whether this version addresses your comment.

[jhenin]

Whenever possible we rebase branches onto master rather than merging master. I've done that here.

[jhenin]

Maybe I was not clear in my first review. Since this is a moving force constant scheme, I think it should integrate into colvarbias_restraint_k_moving and the update can be done in the well-named function update_k().

[jhenin]

What I meant before is, the name of the CV to reference could be provided as a new parameter in the bias configuration, so no reverse lookup would be needed, just a regular lookup by CV name.

I recommend a slightly different separation of tasks between the CV (and its constituent CVC) and the bias.
Since the CV's value is lambda, it does not need to know about the exponent n. The bias needs to know about n, so it should compute the force on lambda for extended Lagrangian integration and apply it to the CV using apply_bias. Effectively, this bias acts on the restrained CVs, but also on the lambda CV, because its energy depends on CV values as well as lambda.

[fhh2626]

Thanks Jerome! So what does your expected Colvars config file look like? Could you please provide an example to prevent misunderstanding? I can revise the code to adapt it.

[jhenin]

The user interface could be something like:

colvar {
  name lambda_cv

  extendedLagrangian on
  extendedMass 1500000
  extendedLangevinDamping 5000
  reflectinglowerboundary  on
  reflectingupperboundary  on

  harmonicForceConstant # Similar to alchLambda, no parameters needed
}

harmonic {
    colvars d
    forceConstant 100.  # k_max
    dynamicForceConstantLambda lambda_cv
    dynamicForceConstantExponent 4.0
}

When initializing the harmonic restraint, it can look up the cv by name and do the linking, that is, set a pointer to itself in the CVC. The CVC would complain if asked to calculate itself when the bias pointer is not set.