Constrained Geometry Optimization in PBC with Gaussian

Fernando R. Clemente, Ph.D.

Technical Support
Gaussian, Inc.

This is from an email exchange between Rich Martin from LANL and Gaussian Technical Support.

The best way to apply constraints in these geometry optimizations would be to use "Opt=ModRedundant", symbolic Z-matrices would not work well for PBC optimizations.

The trick in applying constraints to PBC geometry optimizations is to get the correct numbering right. Note that the translation vectors count for the purpose of atom numbering even though they are obviously not atoms and thus their numbers cannot be used directly to apply a constraint. So, for example, if one wants to freeze the length of a translation vector, the constraint would be to freeze the interatomic distance between an atom in the origin cell and the same atom in the contiguous cell along the direction of the translation vector of interest.

In a 3-D periodic calculation, there are 8 cells to consider, the origin cell (O) and seven contiguous cells (one along each one of the three translation vector, T1, T2 and T3, plus the four contiguous cells along the diagonals, T1+T2, T1+T3, T2+T3, and T1+T2+T3). The following shows how these cells are ordered for atom numbering purposes:

Cell       Start   End
O              1     N
T3           N+1    2N
T2          2N+1    3N
T2+T3       3N+1    4N
T1          4N+1    5N
T1+T3       5N+1    6N
T1+T2       6N+1    7N
T1+T2+T3    7N+1    8N

where N is the number of entries in the input file, that is the number of "real" atoms plus the three translation vectors.

For instance, in the input file you sent here, you have 12 atoms plus the three translation vectors, that makes N=15. If you were to freeze the length of the first translation vector (T1), you could freeze, for instance, the interatomic distance between atom 1 (in cell O) and atom 61 (same atom but in cell T1).

In your case, I see that you would like to constrain the lengths of T1 and T2 to be the same, not just freeze the length. Unfortunately, the redundant internal set does not allow to impose the identity condition between two coordinates. However, for this particular case, it would still be possible to do what you are trying here. In order to keep the lengths of T1 and T2 equal during the optimization, you could freeze the angle between an atom in cell T1, the same atom in cell O, and the same atom in cell T1+T2. This angle should bisect the angle between T1-O-T2. So, given that your translation vectors form an angle of 90 degrees, the T1-O-T1+T2 angle should be frozen to be 45.0 degrees, exactly at the bisection (which implies that the two translation vectors will be equal in length). Below is an example of how to do this for your input file.

#p gen pbc=(nkpoint=2000) pop=regular
scf=(convergence=7,maxcycles=64,novaracc) utpsstpss

V4O8 tetragonal rutile  (high temperature);  Towler V and O basis;

0 5
V                  1.27880629    4.44482660   -0.22374016
O                  0.27860530    0.95052300    0.87756065
O                  1.81959437    3.32683050    0.96076224
V                  3.48744003    2.18167660    1.19205909
V                  2.73205066    0.08147340    3.77512190
V                  0.52341693    2.34462340    2.35932265
O                  4.48764102    3.21367300    0.09075828
O                  3.73225165    3.57577700    2.67382109
O                 -0.47678407    1.31262700    3.46062345
O                  2.94665196    1.06368050    0.00755669
O                  2.19126259    1.19946950    2.59061950
O                  1.06420500    3.46261950    3.54382505
Tv                 4.55460000    0.00000000    0.00000000
Tv                 0.00000000    4.55460000    0.00000000
Tv                 0.00000000    0.00000000    5.70280000

A 61 1 31 90.0 F   ! T1-O-T2 angle 90 deg.
A 61 1 91 45.0 F   ! T1-O-T1+T2 angle bisects above angle


This page was last updated on January 14, 2010 by Cristian V. Diaconu