Relative orientation between two molecules as a CV

[AMPsUSC]

I am trying to get the Standard Free energy for the binding process of two molecules. In my system the relative orientation and the COM distance are the most obvious CVs. I do not see an obvious way of defining the relative orientation. When analyzing MD trajectories of my system I use the dot product between unitary vectors defining the orientation of each molecule. Is there a direct way of doing something similar with the currently implemented CVs?

[mquevill]

Unfortunately, I don't believe that any current CVs in SSAGES can calculate this exact property for you. We are definitely interested in developing CVs based on director vectors. I could imagine an analog of the AngleCV that would calculate angles between molecules' (or groups of atoms) unitary vectors, but we do not have the framework to calculate those types of properties just yet.

[AMPsUSC]

Thank you very much for your answer. Actually I realized I need some other exotic CVs. I thought about using ANNCVs but in my case I know exactly what I want so I would like to consider them via geometrical relationships. I took a look to this presentation http://miccom-center.org/images/MiCCoM_CV_talk_revealed.pdf that explains how to implement new CVs but it seems it was written for a different version of SSAGES since I cannot find the CoolThingCV files in the source code...

[mquevill]

If you have a well-defined property that describes your system, it's probably better to calculate that outright instead of using ANNCVs.

The CoolThingCV that is mentioned in that presentation is a hypothetical CV, and the instructions go through what you need to add to the source tree to properly implement a new CV.

There have been some changes since that presentation was made, but most of the information about adding a CV are the same. You would need to modify src/CVs/CollectiveVariable.cpp, src/CVs/CoolThingCV.h, and schema/CVs/coolthing.cv.json.

I might suggest looking at the changes made in commit https://github.com/SSAGESproject/SSAGES/commit/a156db8c89f922d7b7d1791b79c94c69e81c24b9, since this was the most recently added CV. It only modifies 4 files (including adding documentation), and it shows most of the changes that you would need to add a CV. If it is similar to any existing CVs, I would suggest copying the framework and adapting it to your needs.

If you would like to contribute a new CV that you've added, you could even open a pull request if you think others might benefit from it.

[AMPsUSC]

Opps! :) I got it. Thanks! I will see what I can do.

[hsidky]

i worked on a CV which involved calculating and biasing the orientation of a molecule. This was some years ago so the code is outdated, but I've included the relevant pieces of code below. The key piece is the gradient of the eigenvectors of the inertial tensor. Once you have those you can perform whatever mathematical operation you want on them (like, say, a dot product) then just apply the chain rule. You can see the derivation of the appropriate gradients in the SI of my paper here: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.107801

        std::vector<Vector3> directors;

        for(auto& group : groups_)
        {
            // Get center of mass. 
            auto mtot = snapshot.CachedTotalMass(group);
            Vector3 com = snapshot.CachedCOM(group, mtot);

            // Increment N now, and override later with proper number of selected groups.
            ++N;

            // If not within region, skip.
            if(!IsWithinRegion(com))
                continue;

            selected.push_back(N-1);

            // Initialize inertial tensor. 
            Matrix3 I = Matrix3::Zero();

            // Go through individual local atoms.
            for(auto& id : group)
            {
                auto x = snapshot.ApplyMinimumImage(gpos.row(id).transpose() - com);
                auto& m = mass[id];
                I.noalias() += m*(Matrix3::Identity()*x.squaredNorm() - x*x.transpose());
            }

            // Perform EVD. The columns are the eigenvectors. 
            // SelfAdjoint solver sorts in ascending order. 
            SelfAdjointEigenSolver<Matrix3> solver;
            solver.computeDirect(I);
            const auto& eigvals = solver.eigenvalues();
            const auto& eigvecs = solver.eigenvectors();

            // Assign variables for clarity.
            auto l1 = eigvals[0], 
                 l2 = eigvals[1], 
                 l3 = eigvals[2];
            const auto& u1 = eigvecs.col(0), 
                        u2 = eigvecs.col(1), 
                        u3 = eigvecs.col(2);

            // Make sure the eigenvector is always pointing towards the first atom.
            // Collect first and second atom coordinates. 
                        const Vector3& x0 = gpos.row(group[0]), x1 = gpos.row(group[1]); 
            Vector3 u = u1;
            if(u1.dot(snapshot.ApplyMinimumImage(x0 - x1).normalized()) < 0)
                u = -u1;

            // Chain rule #1. du_i/dx_j
            for(auto& id : group)
            {
                auto x = snapshot.ApplyMinimumImage(gpos.row(id).transpose() - com);
                auto& m = mass[id];
                gradient_[id].noalias() = 
                    1./(l3 - l1)*m*u3*(x.dot(u3)*u + x.dot(u)*u3).transpose() + 
                    1./(l2 - l1)*m*u2*(x.dot(u2)*u + x.dot(u)*u2).transpose();
            }

        }

Note that this code does not work. You will have to re-work it but it should be a starting point.

[AMPsUSC]

Thank you! I had already started to work in a different approach.... but I will see your code later

[AMPsUSC]

One of the CV I want to add is the angle between v1 and v2, being v1 the vector joining the COM of groups 1 and 2 and v2 the vector perpendicular to the plane formed by a third group of atoms. Using my very poor C knowledge, I added a function named PerpendicularUnitaryVector to Snapshot.h to get v2 from a set of atom ID's, while v1 is directly calculated in Evaluate. I hope what I did is correct but now I have doubts on how to implement the calculation of the gradient. My current Evaluate function is as follows:

                void Evaluate(const Snapshot& snapshot) override
                {
                        // Get local atom indices.
                        std::vector<int> idx1, idx2, idx3;
                        snapshot.GetLocalIndices(group1_, &idx1);
                        snapshot.GetLocalIndices(group2_, &idx2);
                        snapshot.GetLocalIndices(group3_, &idx3);

                        // Get data from snapshot.
                        auto n = snapshot.GetNumAtoms();

                        // Initialize gradient.
                        std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
                        grad_.resize(n, Vector3{0,0,0});

                        // Get centers of mass.
                        Vector3 com1 = snapshot.CenterOfMass(idx1);
                        Vector3 com2 = snapshot.CenterOfMass(idx2);
                        Vector3 v2 = snapshot.PerpendicularUnitaryVector(idx3);

                        // Account for pbc.
                        Vector3 v1 = snapshot.ApplyMinimumImage(com2 - com1).cwiseProduct(dim_.cast<double>());
                        v1=v1/v1.norm();

                        dotP=v1.dot(v2);
                        val_ = acos(dotP);
                }

Any help to get the gradient would be appreciated.

[mquevill]

You will have to propagate the gradient through the PerpendicularUnitaryVector, since that calculation will introduce some effects into the gradient. The gradient is with respect to each atom's coordinates.

You can get the gradient of v1 like is done in ParticleSeparationCV. The logic to calculating the gradient from getting the angle between the two vectors can be taken from AngleCV. It can be somewhat tricky to track down where all of the terms in the gradient come from though.

Be sure to propagate gradients through every function. For instance, in CenterOfMass(), you will get a term of mass[i]/masstot due to mass weighting within the COM calculation.

[AMPsUSC]

I got it. A bit tedious in this case but easy to understand. Thank you!

[AMPsUSC]

I implemented the new CV and it seems to be working fine. Thank you for your help.

BTW the atom IDs do not correspond to those of the input pdb file? while debugging it seems that there is a -1 shift in the index list, as if SSAGES starts to count the atoms starting from 0 instead of from 1 (as in standard pdb files).

[mquevill]

Since you were able to get the CV up and running, I will close this issue. We can continue the index discussion in #23.