Sire

The code in Sire to handle complete molecular systems is now progressing well, and I am now finally implementing ideas that I have had way back in April last year. In particular, I have now pretty much finished the implementation of dynamic forcefields. Each forcefield has functions (FFComponents) that represent components within that forcefield, e.g. the coulomb and LJ energy of an intermolecular forcefield, or the bond energy of an intramolecular forcefield. Now I have written the code that allows the total energy of the system to be represented as an arbitrary expression that involves those functions! For example, I can now write;


e_protein_ligand_a = e_coul_protein_ligand_a + e_lj_protein_ligand_a
e_solvent_ligand_a = e_coul_solvent_ligand_a + e_lj_solvent_ligand_a
e_intra_ligand_a = e_total_intra_ligand_a
e_ligand_a = e_protein_ligand_a + e_solvent_ligand_a + e_intra_ligand_a

e_protein_ligand_b = e_coul_protein_ligand_b + e_lj_protein_ligand_b
e_solvent_ligand_b = e_coul_solvent_ligand_b + e_lj_solvent_ligand_b
e_intra_ligand_b = e_total_intra_ligand_b
e_ligand_b = e_protein_ligand_b + e_solvent_ligand_b + e_intra_ligand_b
e_rest_of_system = e_coul_protein_solvent + e_lj_protein_solvent + e_intra_protein

lam = Symbol("lambda")

e_total = e_rest_of_system + (1 - lam)*e_ligand_a + lam*e_ligand_b

lam_f = Symbol("lambda_f")

e_forwards = e_rest_of_system + (1-lam_f)*e_ligand_a + lam_f*e_ligand_b

system.setTotalEnergy( e_total )
system.monitorAverage( e_forwards - e_total, ExponentialAverager() )

In other words, you can completely control how Sire calculates energies (and thus by implication forces and other derivatives of the potential energy), and which energies are averaged. The example above shows how the forcefield for a dual-topology free energy simulation could be implemented (assuming that the user has also linked the motion of ligand A and ligand B together). This shows how Sire itself does not need to know anything about free energy calculations! Thus, pretty uniquely I think, Sire is a free energy code that does not know anything about free energy calculations! The great thing about this is that all of the energy components are calculated rapidly (faster than ProtoMS) and the algebra system I built in is sufficiently clever to recognise that e_rest_of_system cancels out when calculating e_forwards - e_total. I hope that it is clear that this is an incredibally powerful way of building a simulation program, as it allows very rapid prototyping of new ideas, as it is possible to build a forcefield pretty much any way you wish, and that you can monitor the average of pretty much any property you wish (and control the energy using whatever parameters you wish, e.g. lam, lam_f etc.). This could be very powerful....

On an unrelated note - I found out today that Lennard-Jones was the Dean of Science here at Bristol when he presented the work that lead to the Lennard-Jones equation. And the paper was published just over 75 years ago - so just missing an important anniversary. I read the paper today ("Cohesion", J.E. Lennard-Jones, Proceedings of the Cambridge Philosophical Society, Vol. 43, Part 5, Sept. 1st 1931 - google for it or contact me if you want a copy) and it is a very interesting read. Especially so as he doesn't present the actual equation (indeed, the equation is probably something that he would not have expected, as the paper's aim was to use QM to find the distance dependence of attractive and repulsive vdw forces). Amazing.