I'm sure you'll want some context
"but WHYYY are you checking virials in sander??"
So, here you go:
[context]
I want to use sander to move atoms around for me, but I want to compute the energy and forces myself using a different Hamiltonian.
I do this by computing my energy and forces near the top of force.F90 and then skip the evaluation of the "regular" force components in sander, including ewald.
This appears to work fine for minimizations and NVE. I want to do some constant pressure simulations, so I need to compute the virial(s) and then TEST whether or not I computed them correctly.
Before implementing the virials, I want to verify that sander is computing the virials correctly. I mean, regular-sander in git without any modifications by me. It certainly must be computing them correctly, but I need to check them just to verify that my procedure for TESTING is correct or not. This also helps me understand how to insert my virials into sander's data structures.
[/context]
OK.
I'm having difficulties verifying that the virials in regular-ol'-sander are correct.
chkvir in debug.F90 is clearly broken and has been for some time now. The git logs show changes to the force() calls, but the last significant change to the actual testing of virials occurred in ~ Nov 2000, I think. The configuration scripts required to build those old commits don't appear to be in the repository, so I haven't been able to trace back to a point in time when debug.f (at that time) was working.
I've spent the past few days trying to debug debug.F90. You may have noticed that the version in git yields warnings like
"""
debug.F90:326.12:
debug.F90:329.12:
debug.F90:333.12:
debug.F90:336.12:
ene,del,duda,atvir,apduda,type, qsetup)
1
Warning: Type mismatch in argument 'ene' at (1); passed TYPE(state_rec) to REAL(8)
"""
The changes involve the interface to force(), which now use the types defined in the "state" module. The energy is then obtained from ene%pot%tot instead of ene(23)
In addition, the calls in debug_frc() to check_vtens erroneously passed-in molvir and atvir - which are global variables that are overwritten with each call to force() - so I pass in the "saved" versions mvten and avten. Presumably this was the intention when written.
I have attached a patch debug.F90.patch that can be applied to SHA b327556ed17c19a8be629467312b56ce48d1778a
patch -c debug.F90 debug.F90.diff
With the above changes, chkvir will now output nonzero values (hey - that's a start); however, the values are mostly garbage. I mean, the comparison between the analytic and numerical dU/da's are usually wrong.
Although I haven't found a way to make the checks work, I have made a few observations.
Observations:
[These observations are made with a 512 tip3p cube of water. Shake is turned off. I displace the first water's oxygen away from its equilibrium position to guarantee a bond contribution to the energy. I've tried using debugf within NVE and NPT - no difference.]
If my &debugf is
&debugf
do_debugf = 1,
neglgdel = 5
zerochg = 1
zerovdw = 1
do_bond = 1
do_angle = 0
chkvir = 1
/
i.e., the only energy contribution is from the bond energy, then the atomic and molecular virials dU/da tests pass. In this case, the atvir and molvir computed by force() are zero (apparently they are only the nonbond-part of the virial). Therefore, in this case, the debug test compares the numerical dU/da to that computed from the intvir array, which is evaluated in get_intvir():
intvir_{alpha,beta} = -\sum_i f_{i,alpha} r_{i,beta}
If my &debugf is
&debugf
do_debugf = 1,
neglgdel = 5
zerochg = 1
zerovdw = 0
do_bond = 0
do_angle = 0
chkvir = 1
/
i.e., only vdw energy contributes, then the tests fail. They aren't even close. I note, however, that if I change
call check_vtens(xx,ix,ih,ipairs,x,f, &
ene,del,duda,avten,apduda,type, qsetup)
to instead read
call check_vtens(xx,ix,ih,ipairs,x,f, &
ene,del,duda,intvir,apduda,type, qsetup)
(so it is using -\sum_i f_{i,alpha} r_{i,beta}), then all off-diagonal elements of dU/da match finite difference, but the diagonal elements do not.
I'm more interested in trying to get the atomic virials to work first because the molecular virials are related to them.
At this point, I'm kind of stuck and scratching my head as I try to figure out what is going on in the code.
I was hoping to gain a better understanding by seeing how the virials are used in ew_force.F90, but this just raised more questions in my mind. For example, in Essmann, J. Chem. Phys., 103(19), 15 (1995) on pg 8597 (end of top paragraph on right) they state: The molecular virial is the atomic virial MINUS
\sum_{mol} \sum_{i in mol} f_{i,alpha} ( r_{i,beta}-r_{mol_com,beta} )
where f is the total NONBOND force on atom i in direction alpha, but in sander/ew_force.F90 subroutine get_nonbond_virial(), the line
molvir(i,j) = molvir(i,j) + frc(i,n)*xr(j,n)
appears to add in the nonbond forces, instead of subtracting them. Printing out xr and comparing to my input coordinates, it appears that xr really is r_{i,beta}-r_{mol_com,beta}, so it looks like it is the wrong direction from my naive understanding (unless molvir has the opposite sign of what I think it does.)
If anyone has any comments, thoughts, or suggestions, I'd be happy to hear them.. I'm not getting anywhere with what I am doing.
-Tim
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Thu Dec 05 2013 - 17:00:03 PST
