Normal modes with vsites and/or shells does not work
Tested normal mode analysis with different small molecules. TIP5P gives NaN everywhere in the mtx file. TIP4P gives strange numbers as well, while models with shells yield weird numbers. I guess that in all cases mdrun should only work on the real atoms and not on the vsites/shells.
Normal modes don't work currently with virtual sites or shells.
Made normal modes work with shells and vsites
Implement shells and vsites in normal mode analysis (do_nm routine)
and analysis of eigenvalues and frequencies. The normal mode analysis
is done on real atoms only and the shells are minimized at each step
of the analysis.
#1 Updated by Berk Hess over 5 years ago
- Status changed from New to Feedback wanted
This is not trivial to fix. You can not have mdrun only work on the normal atoms.
You need to put in the dependence of the coordinates which are a function of the
mass into the matrix calculation. Herman Berendsen derived a formula for shells.
We need to add a formula for virtual sites, which shouldn't be hard.
Who has time to look into this?
#2 Updated by David van der Spoel over 5 years ago
I will try. Do you have any equations written down? I would say that in both cases the normal mode calculation should be performed on the particles with particle type atom, not vsite or shell. Then, if you displace an atom you have to recompute the virtual sites, and in the case of shells, minimize them. Is it more complex than that? Maybe the .mtx file should only store the atomic hessian, otherwise g_nmeig should be changed too.
#12 Updated by Berk Hess about 2 years ago
I think that for shells the trick is as follows:
You can write the Hessian for the full system of masses+shells as:
where A is the Hessian for the masses and B for the shells. Then the Hessian for the masses coupled by the shells is:
A - C^T B^-1 C
I don't know if this is better than the solution of minimizing the shells for each atom displacement.
#14 Updated by Berk Hess about 2 years ago
For vsites you would need to figure out matrix C for every construction type, which is tedious. Since the vsite construction is exact (and I assume handled correctly) in your current change, there is no reason to use the matrix trick. For vsites there is the iterative minimization, which requires setting and extremely low tolerance, as low as the tolerance used to minimize the structure, and it requires iterations that could increase the cost of the matrix calculation by an order of magnitude. But I guess the matrix calculation is still cheap. The matrix trick would avoid those two issues, but both are not really problematic.
#15 Updated by David van der Spoel about 2 years ago
OK, now that this passed the build tests in Jenkins we should consider whether we want to put much more effort in this. I can build in a warning about the tolerance for shells. I have some test cases but do not know how easy that is to add to the regression tests seeing that you can only run this really in double precision.