Computing normal modes in angle-axis
On this page I will describe how you can compute normal modes in angle-axis framework. It is likely that you may have to implement your own potential (until rigid body sandbox works!), so I will cover that too.
I presented the maths behind normal mode calculations in one of our group meetings. The slides are here: talk
Here are the OPTIM fortran files that are required in my normal modes implementation. In all cases, search for "hk286" and/or "RBAAT" for the relevant changes.
1. modhess.f90 - added a variable which toggles between moving and stationary frames
2. keywords.f - additional initialisation
3. geopt.f - subroutine GEOPT calls VDUMP, which dumps the normal modes and their corresponding eigenvectors. Relevant OPTIM keywords are DUMPVECTORS and ENDHESS.
4.path.f & CONNECT/ncutils.f90 - both files seem to be doing very similar things. Probably different implementation of connect routines. For normal modes calculation purposes, they write to path.info the eigenvalues at the minima and transition states in triplets. Relevant OPTIM keyword is NOFRQS.
5. rigidb.f90 - subroutine NRMLMD is what computes the normal modes. Additionally subroutine COMPUTEINERTIA is a wrapper that points to the right subroutine for computing the inertia matrix of your rigid body / molecule.
For an example that is working look at newtip.f90.
If you need to add a new potential, here are the steps you need to do:
1. Write a subroutine for computing the energy, gradient and hessian. Look at subroutine NEWTIP in newtip.f90 as an example. In particular, make sure you allow switching between moving and stationary frames. This is done by looking up the value of RBAANORMALMODET.
2. Define your rigid body / molecule. See DEFTIP4 in newtip.f90 as an example.
3. Write a subroutine to compute the inertia tensor for the molecule given a rotation matrix RMI. You also need to output the total mass of your molecule. See INERTIANTIP in newtip.f90 as an example. Edit subroutine COMPUTEINERTIA to add the call to your new inertia routine.
4. If you are feeling brave, you can write a visualisation routine to see your normal modes. The math is rather hairy. So look up VISUALISEMODESTIP4 in newtip.f90 as an example. In that block of code, I explained the steps. At the moment, this is still in development stage. So you need to wire it in somewhere if you want it. I suggest putting it on vdump.f.
Getting the numbers out from OPTIM runs:
1. Use ENDHESS
2. Use DUMPVECTORS and DUMPMODE - the output is vector.dump. It gives you the normal mode frequencies followed by the corresponding eigenvectors in angle-axis coordinates.
3. DO NOT USE NOFRQS if you want the normal modes frequencies are printed in path.info for the minima and transition states.
What do the numbers correspond to?
In my implementation, the printed numbers are the frequencies, NOT frequencies squared. Negative frequency actually means imaginary frequency. Unless it is very close to zero, an imaginary frequency implies you are not at a minimum. Unit of the frequency: cm^{-1}.