Optimizing (nonstationary) conformations along a manual dihedral angle scan

Hello,

I am having a difficult time optimizing (via HF/6-31G(d) with g_convergence gau_tight) conformations with fixed dihedrals. As a test system, I am using the heavy atom rotation of butane. The local maxima and minima along the internal coordinate do optimize (i.e. CCCC=0, 60, 120 & 180 degrees). However, the intermediate states do not converge (i.e. CCCC=30, 90 & 150). I have tried different combinations of including/omitting the following input flags: “set intrafrag_step_limit 0.1”, “dynamic_level = 1”, and "set step_type nr ". I have also tried to use the cc-pVTZ basis set, but to no avail. Any suggestions on how I might obtain convergence of these conformations? Below is an example with butanes’s CCCC torsion at 30 degrees.

Psi4 version = Psi4 1.1a2.dev170

Cheers,
Karl

P.S. On a side note, WebMO Version 17.0.008 can not read the file output of Psi4 1.1a2.dev170 (it could read Psi4 1.0.0) of successfully optimized files. It reports “Invalid File Format.”


memory 30 Gb

molecule molecule_2 {
0 1
C 1.553043914483 -0.380563070329 0.177606221727
H 2.315033007414 -1.152987247068 0.193492310224
H 1.732054058031 0.238473154563 -0.693962571225
H 1.697080893325 0.239870251822 1.058865080134
C 0.159486519325 -1.010938954953 0.137512577354
H 0.051990622598 -1.722212020210 0.951236813969
H 0.094522780862 -1.595565747029 -0.778434460525
C -1.020534480320 -0.018198532697 0.172889188287
H -1.886554491249 -0.478069985632 -0.294071439145
H -1.301033945721 0.156183601950 1.210195327922
C -0.746887241731 1.344573818388 -0.467118028830
H -1.641488698365 1.958522401339 -0.444992962509
H 0.034938821002 1.883681666067 0.055248446015
H -0.442961882564 1.247556855232 -1.506309702384
}
set globals {
reference rhf
g_convergence gau_tight
basis 6-31G(d)
intrafrag_step_limit 0.1
}

set optking {
fixed_dihedral = ("
1 5 8 11 30.0
")
dynamic_level = 1
}

optimize(‘scf’)

Of course beginning from the best starting structure you can is helpful. Try increasing the value of FIXED_COORD_FORCE_CONSTANT (its default is 0.5). The progressively increased force constant applied around the desired torsional value is determined by the formula:

k = (1 + 0.05 * (iter-1)) * fixed_coord_force_constant

so you could try a substantially larger value if you are close to the target but not quite getting there. This formula is a compromise for cases where the initial steps need to be gentle, and is not ideally suited for very tight convergence criteria.

Thanks Rollin for the feedback. By setting FIXED_COORD_FORCE_CONSTANT to 160.0 I was able to optimize butane’s constrained 30, 90 and 150 degree angle (CCCC) using the gau_tight convergence criteria. Note that this value was found by iteratively increasing it from 0.5 to 160.0 (e.g. 0.5, 1.0, 1.5, 3.0, 10.0, etc). Consequently, the specification of intrafrag_step_limit and dynamic_level in my example above was not needed.

I can also confirm that the resulting relative energy profile and geometries are in very good agreement with the results that Gamess provides. (Small note: the Gamess calculations were done without using DF.)

Bests,
Karl

Hi Karl, how does the working input look after your corrections? I am trying to run this command for a dicyclic molecule, and make it rotate about the center. How do I define the center bond to rotate about?

Thanks

Hello Sergio,

Below is the final version of the input that works for me. One can add the other keywords if the optimization gives you trouble: “intrafrag_step_limit 0.1” and “dynamic_level = 1”. You define the torsion angle using fixed_dihedral = ("1 5 8 11 0.0 ") where the first four integers are the sequential atom numbers that define the angle of interest (i.e. C1-C5-C8-C11), and the last number is the angle value (i.e. zero degrees).

Bests,
Karl

memory 30 Gb

molecule molecule_1 {
0 1
C          1.95400        0.13200       -0.00000
H          2.73800       -0.61900        0.00000
H          2.09600        0.75700       -0.87700
H          2.09600        0.75800        0.87700
C          0.56800       -0.51200        0.00000
H          0.47000       -1.15700        0.87100
H          0.47000       -1.15800       -0.87000
C         -0.56800        0.51200       -0.00000
H         -1.20000        0.34800       -0.87000
H         -1.19900        0.34800        0.87100
C         -0.07100        1.95700        0.00000
H         -0.89900        2.65900        0.00000
H          0.53700        2.16300        0.87700
H          0.53600        2.16300       -0.87700
}
set globals {
    reference rhf
    g_convergence gau_tight
    basis 6-31G(d)
    fixed_coord_force_constant 160.0
}
 
set optking {
    fixed_dihedral = ("
    1 5 8 11 0.0
    ")
}

optimize(‘scf’)

Thanks Karl. I will try your new file on the system here.

Cheers

Hi Karl, maybe this error is not related to the command itself, but I used avogadro to find the atom numbers in the dihedral, and then inserted them as you did, and here is what I got:

memory 20 Gb
molecule {
0 1
  Cl       -0.93360       -1.72050       -1.77550
  Cl        0.93360       -1.72040        1.77560
  Cl       -0.22090        2.13990        2.01800
  Cl        0.22100        2.13970       -2.01810
  Cl       -4.06040       -1.66870       -1.11090
  Cl        4.06040       -1.66870        1.11090
   C       -0.71460        0.20170        0.13980
   C        0.71460        0.20170       -0.13980
   C       -1.57020       -0.64930       -0.55970
   C        1.57010       -0.64920        0.55980
   C       -1.22650        1.05920        1.11360
   C        1.22650        1.05910       -1.11370
   C       -2.93780       -0.64270       -0.28550
   C        2.93770       -0.64270        0.28560a
   C       -2.59420        1.06570        1.38780
   C        2.59420        1.06560       -1.38780
   C       -3.44980        0.21490        0.68820
   C        3.44980        0.21480       -0.68820
   H       -3.01030        1.72740        2.14330
   H        3.01030        1.72730       -2.14350
   H       -4.51280        0.23360        0.91520
   H        4.51290        0.23350       -0.91530
}
set globals {
reference rhf
g_convergence gau_tight
basis cc-pVTZ
fixed_coord_force_constant 160.0
}

set optking {
fixed_dihedral = ("
1 2 4 5 180.0")
}
optimize('CCSD(T)')



--------------------------------------------------------------------------

Traceback (most recent call last):
  File "/usr/bin/psi4", line 248, in <module>
    exec(content)

SyntaxError: invalid syntax (<string>, line 52)

*** Psi4 encountered an error. Buy a developer more coffee!
*** Resources and help at github.com/psi4/psi4.

Hi again,

Part of the problem with how I was getting the text formatted for viewing on the webpage. I have corrected it above. Note that indentations are important in python.

However, your input also include an “a” at the end of a carbon coordinate line. You also need to move the the ") after 180.0 to another line that is indented (see above).

I would also recommend defining the torsion between consecutively linked atoms. For your goal, I believe atoms 9 7 8 10 would do. You may also want to seriously reconsider doing an optimization at the CCSD(T)/cc-VTZ level - this is very very large for the number of electrons your system has.

Good luck.

Bests,
Karl

Speaking of syntax, I edited @sem75’s last post so the Psi file is enclosed in triple backticks (`). That makes it play nicely with the forums.

Two other things to point out:

  1. In Psi 1.3, the program should tell you what the syntax error is when you have one. Let us know if that isn’t the case.
  2. As Karl said, CCSD(T)/cc-VTZ would be extremely expensive for your system, and especially in Psi4. First, Psi4 doesn’t have analytic gradients with the frozen core approximation. (If you didn’t know, all electron computations are the default in Psi.) Second, CCSD(T) gradients in Psi are, sadly, much slower than in other codes.