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

[Axel Kohlmeyer]

On Mon, 7 Jan 2008, Philip Shemella wrote:

DK> Dear CPMD'ers,

dear daria (& phil and the rest),

not sure whether anybody cared to respond to this so far.
i had the mail sitting in my CPMD inbox as a potential
"look into it on occasion" issue for quite a while.

i did not yet look into all problems listed below, but 
for some other reasons i just had a look into how constraint
and restraint forces and energies are computed and collected,
... and found a number of inconsistencies. depending on the
type of calculation that you do, restraints may actually be
added _twice_ and the combined constraint/restraint energy may 
not be properly reset for all kinds of calculations. worst
case would be to have restraints with no constraints in a
geometry optimization or born-oppenheimer MD.

i'm currently trying to find a solution that will be consistent 
across all types of calculations, but due to the nature of
the affected parts of the CPMD code it is very likely, that 
there is no simple fix, i.e. anybody affected might be forced
to wait until the next release or may need to add specific 
"ad-hoc modifications" depending on the job type which will
break the code for other job types. since this would be ugly 
hacks, i don't want to make them public. please contact me 
directly, if you see some oddities with contraints/restraints 
for what changes you may try to get the correct behavior.

cheers,
    axel.

DK> 
DK> I'm relaying this message for a friend who is not yet on the mailing 
DK> list -- please reply to Daria at dkhvosti at uiuc.edu
DK> 
DK> Best,
DK> Phil Shemella
DK> RPI
DK> 
DK> --
DK> 
DK> Dear CPMD subscribers,
DK> 
DK> I am trying to understand how to run structure
DK> optimizations with restraints. I have read the
DK> manual, but I still have several questions
DK> regarding the units of force constants and the
DK> values of forces resulting from such
DK> calculations, particularly for angle bends.
DK> 
DK> I have run several gradient calculations on
DK> optimized methane molecule applying an angle bend
DK> restraint for one of the H-C-H angles. Molecular
DK> geometry (coordinates of atoms) was the same in
DK> all calculations. Gradients were zero (below
DK> threshold for geometry optimization) in absence
DK> of restraints since geometry was optimized.
DK> 
DK> In the tests the value of the force constant kval
DK> was the same in all calculations, but the value
DK> of the equilibrium angle theta_eq varied.
DK> 
DK> My understanding is that in a system like this
DK> the strain energy of the restraint is
DK> E=0.5*kval(theta –- theta_eq)^2, and forces on
DK> hydrogen atoms participating in the restraint are
DK> related to the deviation from the equilibrium angle as follows:
DK> Force = kval*(theta -– theta_eq)/distance(H-C).
DK> (By Force I mean the magnitude of the force vector).
DK> 
DK> Based on that relationship, the value of Force
DK> should be proportional to the deviation of the
DK> bending angle from the equilibrium angle if other
DK> variables are unchanged. However, the plot of
DK> Force vs. (theta –- theta_eq) showed quadratic dependence.
DK> 
DK> I would appreciate if somebody could clarify why
DK> magnitude of the force vector is not directly
DK> proportional to (theta - theta_eq) in these
DK> calculations. Also, I am not sure what the units
DK> of kval in the input are. Is it “Hartree/rad^2”?
DK> 
DK> Just in case, the relevant portion of the input
DK> file is below. I used BLYP functional and input
DK> in atomic units in all calculations; values of
DK> theta_eq were 30-80 (in increments of 10); I also ran a test with 
DK> theta_eq=0.
DK> 
DK> &SYSTEM
DK>    SYMMETRY
DK>      1
DK>    CELL
DK>          27.8182702875 1.0 1.0 0.0 0.0 0.0
DK>    CUTOFF
DK>      30.0
DK> &END
DK> 
DK> &ATOMS
DK> *C_VDB_BLYP.psp NEWF FORMATTED
DK>     LMAX=P
DK>       1
DK>    -2.3418595788 2.7361917529    -0.0861679775
DK> *H_VDB_BLYP.psp NEWF FORMATTED
DK>     LMAX=S
DK>       4
DK>          -2.4048529099   0.6661909468    -0.1372139660
DK>          -1.3569794985   3.4388324334    -1.7678168469
DK>          -4.2730035025   3.4853022511    -0.0536468133
DK>          -1.3325198437   3.3546522665    1.6139082542
DK> 
DK>   CONSTRAINTS
DK>    RESTRAINTS
DK>     1
DK>   BEND 2 1 3 theta_eq 5.000
DK>   END CONSTRAINTS
DK> 
DK> &END
DK> 
DK> 
DK> Thank you very much.
DK> 
DK> 
DK> _______________________________________________
DK> CPMD-list mailing list
DK> CPMD-list at cpmd.org
DK> http://cpmd.org/mailman/listinfo/cpmd-list
DK> 

-- 
=======================================================================
Axel Kohlmeyer   akohlmey at cmm.chem.upenn.edu   http://www.cmm.upenn.edu
   Center for Molecular Modeling   --   University of Pennsylvania
Department of Chemistry, 231 S.34th Street, Philadelphia, PA 19104-6323
tel: 1-215-898-1582,  fax: 1-215-573-6233,  office-tel: 1-215-898-5425
=======================================================================
If you make something idiot-proof, the universe creates a better idiot.