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

[Philip Shemella]

Dear CPMD'ers,

I'm relaying this message for a friend who is not yet on the mailing 
list -- please reply to Daria at dkhvosti at uiuc.edu

Best,
Phil Shemella
RPI

--

Dear CPMD subscribers,

I am trying to understand how to run structure
optimizations with restraints. I have read the
manual, but I still have several questions
regarding the units of force constants and the
values of forces resulting from such
calculations, particularly for angle bends.

I have run several gradient calculations on
optimized methane molecule applying an angle bend
restraint for one of the H-C-H angles. Molecular
geometry (coordinates of atoms) was the same in
all calculations. Gradients were zero (below
threshold for geometry optimization) in absence
of restraints since geometry was optimized.

In the tests the value of the force constant kval
was the same in all calculations, but the value
of the equilibrium angle theta_eq varied.

My understanding is that in a system like this
the strain energy of the restraint is
E=0.5*kval(theta –- theta_eq)^2, and forces on
hydrogen atoms participating in the restraint are
related to the deviation from the equilibrium angle as follows:
Force = kval*(theta -– theta_eq)/distance(H-C).
(By Force I mean the magnitude of the force vector).

Based on that relationship, the value of Force
should be proportional to the deviation of the
bending angle from the equilibrium angle if other
variables are unchanged. However, the plot of
Force vs. (theta –- theta_eq) showed quadratic dependence.

I would appreciate if somebody could clarify why
magnitude of the force vector is not directly
proportional to (theta - theta_eq) in these
calculations. Also, I am not sure what the units
of kval in the input are. Is it “Hartree/rad^2”?

Just in case, the relevant portion of the input
file is below. I used BLYP functional and input
in atomic units in all calculations; values of
theta_eq were 30-80 (in increments of 10); I also ran a test with 
theta_eq=0.

&SYSTEM
   SYMMETRY
     1
   CELL
         27.8182702875 1.0 1.0 0.0 0.0 0.0
   CUTOFF
     30.0
&END

&ATOMS
*C_VDB_BLYP.psp NEWF FORMATTED
    LMAX=P
      1
   -2.3418595788 2.7361917529    -0.0861679775
*H_VDB_BLYP.psp NEWF FORMATTED
    LMAX=S
      4
         -2.4048529099   0.6661909468    -0.1372139660
         -1.3569794985   3.4388324334    -1.7678168469
         -4.2730035025   3.4853022511    -0.0536468133
         -1.3325198437   3.3546522665    1.6139082542

  CONSTRAINTS
   RESTRAINTS
    1
  BEND 2 1 3 theta_eq 5.000
  END CONSTRAINTS

&END


Thank you very much.