Inconsistency in crosstacking energy

[cabb99]

As a sanity check I created a simulation using the reverse complement of the dna sequence in the bdna_curv example. Then I reran the example simulation using the reverse complement system hamiltonian. I took care that every atom was in the position expected on the reverse complement system (Atoms in chain B should be before the atoms in the chain A). If the hamiltonian is consistent, the energies should be exactly the same for both cases. During this experiment most of the energies were exactly the same, with the exception of the crosstacking energy.

The base pair energy in the reverse complement hamiltonian is similar to the base pair energy in the Original Hamiltonian. Other energies, such as bond, angle, dihedral, exclusion and electrostatics are also consistent [not shown].

The crosstacking energy in the reverse complement hamiltonian is different to the crosstacking energy in the Original Hamiltonian

I think the reason for this difference is that when the sense chain is changed to the reverse complement, (for example setting the sense chain as B and the anti-sense chain as A), the definition for the crosstacking angle θCS also changes. The crosstacking angles definition. θCS is defined using the particles a b e, and c d f respectively. Since b<d usually b is in the sense chain and d is in the anti-sense chain. The crosstacking angles definition in the reverse Hamiltonian. θCS will be now defined using the particles a b g, and c d h respectively. Since the angle θCS is now defined with different particles we can expect θCS to have a different value.

I believe the energies obtained from the 3SPN2.C Hamiltonian shouldn't depend on the order on which the chains are loaded, so in my opinion this is a issue that should be fixed. I attach the simulations I ran on this file. They also contain a script used to transform the trajectory from the Original system, to the Reverse Complement Hamiltonian.

The part of the code that enforces that b should be less than d is in Pair3spn2::compute

                    if (steb < sted)
                    {
                        stee = sted - 3;
                        stef = steb + 3;
                        stea = steb - 1;
                        stec = sted - 1;
                        myitype = itype;
                        myjtype = jtype;
                    }
                    else
                    {
                        int tmp = sted;
                        sted = steb;
                        steb = tmp;
                        stea = steb - 1;
                        stec = sted - 1;
                        stee = sted - 3;
                        stef = steb + 3;
                        myitype = jtype;
                        myjtype = itype;
                    }

[lequieu]

Thanks for looking into this. Please forgive the delay in response, I haven't run 3SPN2 in several years and I had to dust off the cobwebs.

I've thought a bit more about this and I'm still not sure if this is a "bug" that should be fixed or peculiarity of the model. Would you be willing to run a few more tests on different DNA sequences (specifically the end base-pairs)? I'm curious if its these end bases that are causing this discrepancy when you flip the strands.

Depending on what you find, perhaps there's a better solution than your pull request #5