[CPMD-list] force values and constants with constraint & restraint

[Daria Khvostichenko]

Dear Axel,

Thank you for your reply. I actually looked into the source code, and 
found most of the answers I was looking for. Apparently, the 
potential energy functions for various distortions were not defined 
exactly in the way I thought they were, and, consequently, the 
behavior of derivatives (forces) seemed erratic to me. For example, 
the stretch function is defined as E = 0.5 k( r**2 -r_eq**2)**2, 
whereas I thought it was r = 0.5 k ( r - r_eq)^2, which is instead 
the definition of the distant restraint in the code. The distant 
restraint, however, did not work in my tests, neither with geometry 
optimizations, nor with CPMD simulations: in the output file it 
appears that it only picks the first atom out of the two that define 
the restraint. The angle bend restraint is apparently defined as E = 
0.5 k (cos(theta) - cos(theta_eq))**2, i.e. through cosines, rather 
than simply through angles. It would be great if the manual listed 
equations for restraints rather than text descriptions. Anyway, they 
were relatively easy to locate in the source code, which, thankfully, 
did have these equations explicitly written in the comments.

Forces produced by the code in gradient (single-point) calculations 
with restraints applied to a previously optimized molecule were 
identical in magnitude and direction to those calculated manually 
from derivatives of these potential energy functions with k values in 
Hartree/(Bohr^4) for stretch and in "Hartree/rad^2" for angle bend 
(however, the equilibrium angle theta_eq should be in degrees in the 
input file).

Thanks for the warning about the improper behavior of restraints in 
certain types of calculations. I've also found earlier posts about 
double-counting of restraints in BOMD calculations, and about missing 
conversion of units for some restraints.

As for collection of restraint energies, at least in geometry 
optimizations those energies do not seem to be added to the total 
energy of the system, and there doesn't seem to be a separate entry 
for the restraint energy in the output file.

You mentioned that the worst case for restraints is geometry 
optimization with restraints in absence of constraints. Why is that?

Regards,

Daria