SireAdvanced Molecular Simulation Framework

Now that the core code is complete, I've been adding in some useful functionality. Two recent additions are soft-core potentials and constraints.

Soft-core provides a way to soften intermolecular interactions (predominantly coulomb and lennard jones) so that sampling can be improved when atoms are "switched off", e.g. as occurs during a dual topology free energy simulation. Implementing soft-core in ProtoMS (my old code) was a little difficult (according to Julien - who did all the work) and meant hacking the code all over the place. Implementing soft-core in Sire took an afternoon, and only required the addition of a SoftCLJPotential (and InterSoftCLJPotential). These were then functionalised using the template Inter2B3DFF and InterGroup2D3DFF classes to automatically create the InterSoftCLJFF and InterGroupSoftCLJFF forcefields. And it all worked!

Then, as the alpha parameter of the soft-core forcefield must be mapped to the lambda scaling parameter, I've added constraint classes that allow you to add constraints to the System. for example, here's the (useful) code used to create everything needed for a dual topology soft-core simulation;

    # Create the solute-solvent MM forcefield (intermolecular solute-solvent energy)
    mmff0 = InterGroupSoftCLJFF("solute0-solvent MM")
    mmff0.setCoulombPower(6)
    mmff0.setLJPower(3)

    mmff1 = InterGroupSoftCLJFF("solute1-solvent MM")
    mmff1.setCoulombPower(6)
    mmff1.setLJPower(3)

    lam = Symbol("lambda")

    e_mm0 = solventff.components().total() + mmff0.components().total() + \
            solute_intraff.components().total() + solute_intranb.components().total()

    e_mm1 = solventff.components().total() + mmff1.components().total() + \
            solute_intraff.components().total() + solute_intranb.components().total()

    e_total = ((1 - lam) * e_mm0) + ( lam * e_mm1 )

    system.add(mmff0)
    system.add(mmff1)

     # Add constraints so that alpha = lambda for solute 0
    # and alpha = 1 - lambda for solute 1
    system.add( PropertyConstraint( "alpha", FFName("solute0-solvent MM"), lam ) )
    system.add( PropertyConstraint( "alpha", FFName("solute1-solvent MM"), 1 - lam ) )
 

Note how the dual-topology Hamiltonian is specified in full in the python script, as is also the relationship between lambda and the alpha soft-core scaling parameter.

Cool :-)

*UPDATE* - here is the link to a Sire python script that performs a soft-core dual topology free energy calculation between ethane and methanol in a box of TIP4P water.