Scan problem with symmetry

Hi y’all, I am trying to do a BH3NH3 potential energy surface scan using a zmat of fixed symmetry, and I am running into problems. Here is the input file with no symmetry:

molecule bh3nh3 {
0 1
N
B 1 R
H 1 hn3 2 hnb3
H 2 hb4 1 hbn4 3 dih4
H 1 hn5 2 hnb5 4 dih5
H 2 hb6 1 hbn6 3 dih6
H 1 hn7 2 hnb7 4 dih7
H 2 hb8 1 hbn8 3 dih8

hn3 = 1.004382
hnb3 = 110.832
hb4 = 1.217824
hbn4 = 104.478
dih4 = -59.966
hn5 = 1.004489
hnb5 = 110.744
dih5 = 59.986
hb6 = 1.217435
hbn6 = 104.504
dih6 = -179.955
hn7 = 1.004416
hnb7 = 110.900
dih7 = 179.970
hb8 = 1.217541
hbn8 = 104.533
dih8 = 60.022
}

set {
basis cc-pvdz
}

Rval = [1.5, 2.0, 2.5, 3.0, 3.5]
PES =

for d in Rval:
bh3nh3.R = d
r_idx = "1 2 "
set optking frozen_distance = $r_idx
E = optimize(‘scf’)
PES.append((d, E))

print “\n\tcc-pVDZ SCF energy as a function of R\n”
for point in PES:
print “\t%5.1f%20.10f” % (point[0], point[1])

And this is for the one with fixed symmetry:

molecule bh3nh3 {
0 1
B
N 1 R
X 1 1.0 2 90.0
H 2 Rn 1 An 3 dih4
H 2 Rn 1 An 3 dih5
H 2 Rn 1 An 3 dih6
H 1 Rb 2 Ab 3 dih7
H 1 Rb 2 Ab 3 dih8
H 1 Rb 2 Ab 3 dih9

R = 1.5002
Rn = 1.217541
Rb = 1.004446
An = 104.501526
Ab = 110.819753
dih4 = 120.000
dih5 = 240.000
dih6 = 0.000
dih7 = 180.000
dih8 = 300.000
dih9 = 60.000
}

set {
basis cc-pvdz
}

Rval = [1.5, 2.0, 2.5, 3.0, 3.5]
PES =

for d in Rval:
bh3nh3.R = d
r_idx = "1 2 "
set optking frozen_distance = $r_idx
E = optimize(‘scf’)
PES.append((d, E))

print “\n\tcc-pVDZ SCF energy as a function of R\n”
for point in PES:
print “\t%5.1f%20.10f” % (point[0], point[1])

Which is essentially the same as the one before, with a symmetric zmat input.

The problem is that in the symmetric one, the distance R remains as 1.5 A for all the optimizations.
Even though the output explicitly states:

Setting geometry variable R to 2.000000

A noticed a difference between the two output files. The symmetric one that repeats the optimization
for R = 1.5 prints the optimized geom in cartesian coordinates while the one without symmetry prints it
in z-mat format.

  Writing optimization data to binary file.
    Final energy is    -82.6156165478973
    Final (previous) structure:
    Cartesian Geometry (in Angstrom)
        N    -0.0000030307   0.0000001213  -0.8886031343
        B    -0.0000030307   0.0000001213   1.1113968657
        H    -0.0000030307   0.9405698114  -1.2466722075
        H     1.0324255989   0.5960925854   1.3152637417
        H     0.8145938560  -0.4702636032  -1.2466099041
        H     0.0000213569  -1.1921509487   1.3152872650
        H    -0.8145224186  -0.4703089845  -1.2467263196
        H    -1.0324401473   0.5960581299   1.3153145147
    Saving final (previous) structure.
                    --------------------------
                     OPTKING Finished Execution
                    --------------------------

Final optimized geometry and variables:
Molecular point group: c1
Full point group: C1

Geometry (in Angstrom), charge = 0, multiplicity = 1:

N
B             1           R
H             1         hn3      2        hnb3
H             2         hb4      1        hbn4      3        dih4
H             1         hn5      2        hnb5      4        dih5
H             2         hb6      1        hbn6      3        dih6
H             1         hn7      2        hnb7      4        dih7
H             2         hb8      1        hbn8      3        dih8

R         =    2.0000000000
dih4      =  -59.9991876724
dih5      =   59.9984907880
dih6      = -179.9988279123
dih7      =  179.9983793358
dih8      =   60.0008257225
hb4       =    1.2094613694
hb6       =    1.2094608176
hb8       =    1.2094601925
hbn4      =   99.7041082120
hbn6      =   99.7052432317
hbn8      =   99.7065579430
hn3       =    1.0064217820
hn5       =    1.0064218328
hn7       =    1.0064217327
hnb3      =  110.8415202118
hnb5      =  110.8377238811
hnb7      =  110.8448176060

Cleaning optimization helper files.
Setting geometry variable R to 2.500000
gradient() will perform analytic gradient computation.

And this is the output for the computation with the symmetric z-mat:

    Writing optimization data to binary file.
    Final energy is    -82.6158751587106
    Final (previous) structure:
    Cartesian Geometry (in Angstrom)
        B     0.0000000532   0.8238656988   0.0000000000
        N    -0.0000000436  -0.6761343012   0.0000000000
        H     0.4674733216  -1.0541925521  -0.8096876617
        H     0.4674733216  -1.0541925521   0.8096876617
        H    -0.9349468486  -1.0541924624   0.0000000000
        H     1.1705866699   1.1857528375   0.0000000000
        H    -0.5852932200   1.1857529507  -1.0137577271
        H    -0.5852932200   1.1857529507   1.0137577271
    Saving final (previous) structure.
                    --------------------------
                     OPTKING Finished Execution
                    --------------------------

Final optimized geometry and variables:
Molecular point group: cs
Full point group: C1

Geometry (in Angstrom), charge = 0, multiplicity = 1:

B            0.000000053217     0.823865698805     0.000000000000
N           -0.000000043604    -0.676134301195     0.000000000000
H            0.467473321634    -1.054192552071    -0.809687661705
H            0.467473321634    -1.054192552071     0.809687661705
H           -0.934946848585    -1.054192462415     0.000000000000
H            1.170586669865     1.185752837539     0.000000000000
H           -0.585293220018     1.185752950665    -1.013757727063
H           -0.585293220018     1.185752950665     1.013757727063

R         =    2.0000000000

    Cleaning optimization helper files.

Setting geometry variable R to 2.500000
gradient() will perform analytic gradient computation.

The energy however is the same for all the points

cc-pVDZ SCF energy as a function of R

      1.5      -82.6158751587
      2.0      -82.6158751587
      2.5      -82.6158751587
      3.0      -82.6158751587
      3.5      -82.6158751587

Meaning that the R distance is not increasing during the for loop. I don’t know what is
going on so any help is greatly appreciated.

Hi there! I am currently having a similar problem. Did you ever figure this out?

There are a lot of deeper issues with fixed/frozen internal coordinate syntax and documentation at play involved. I’ll defer to @Rollin_King to talk about the future of those keywords.

For an immediate fix, you can use this code as an example:

CN_distances = [1.460 + i * 0.01 for i in range(200)]

for r in CN_distances:
    r_string = "1 9 " + str(r)
    set optking fixed_distance = $r_string
    optimize('b3lyp')

If that doesn’t solve your problem, please post an output file, enclosed in backticks.