Convergence errors during dihedral scan

Geometry Optimization
1

Hello. I’m a molecular dynamics researcher trying to do a dihedral scan on a cyclohexyl acetate dihedral to eventually parameterise in MD (that is to say my QM knowledge is quite lacking).

I seem to have successfully managed to begin rotating and optimising the geometry at fixed intervals of 10 degrees. Seeing that the dihedral (D(2,4,5,10)) had an angle of 153 in the initially optimised geometry, I started at 150 and began adding 10 each cycle, running:

molecule = psi4.geometry("""
24
CHA.3.xyz
 C   -3.443370463455    0.647769656644    0.523542950280
 C   -1.959792017625    0.803964008843    0.703160331616
 O   -1.440827881011    1.679894908141    1.314167916183
 O   -1.296886858140   -0.168443235859    0.090654887976
 C    0.132792874268   -0.186700385789    0.167190539733
 C    0.737221850414    0.717431911434   -0.900072828218
 C    2.263890848138    0.601938019365   -0.910288899292
 C    2.718587007872   -0.849309082851   -1.078870699323
 C    2.089559565477   -1.755022848729   -0.017683497460
 C    0.563692287603   -1.635031555035   -0.008107529092
 H   -3.960379432774    1.451157007695    1.040215076678
 H   -3.757967278367   -0.319461737706    0.916057660454
 H   -3.683766745564    0.666052536222   -0.539809313296
 H    0.422313924444    0.172445636592    1.153883993526
 H    0.432840937126    1.749073688730   -0.720655923361
 H    0.334130586384    0.419982319792   -1.873857459669
 H    2.662872021991    1.000338867476    0.029278427713
 H    2.675192114576    1.224756843754   -1.708071556679
 H    3.808639452287   -0.911760818497   -1.028868558274
 H    2.430318344225   -1.205356867369   -2.074522820986
 H    2.481343810016   -1.481953560380    0.968748324716
 H    2.376730678543   -2.795566543583   -0.186401470489
 H    0.132170590827   -2.241898250449    0.790828821700
 H    0.153168828402   -2.010315508120   -0.950552627170

""")


set {
    basis cc-pvdz
    scf_type df
    intrafrag_step_limit 0.25
    g_convergence gau
}

def define_angles(psi=0):
    angles = np.arange(0, 360, 10) - 170
    idx = np.where(angles == int(psi))[0][0]
    return(np.concatenate([angles[idx:], angles[:idx]]))

# C-dihedral has initial angle of 153, so we start with 150...
Cdih_range = define_angles(150)
PES_C = []
for psi in Cdih_range:
    dih_str = f"2 4 5 10 {psi} {psi}"
    set optking {
        ranged_dihedral = $dih_str
        max_force_g_convergence = 3.0e-3 # x10 from default
        rms_force_g_convergence = 1.5e-03 # x5 from default
    }
    E = optimize('scf')
    out = f"CHA.C.{psi}.xyz"
    print(psi, E)
    PES_C.append((psi, E))
    molecule.save_xyz_file(out,1)

psi4_dih_scan.dat (2.1 KB)
(I upped the default force values just to see if I could get the code to run through the scan first - rookie mistake?)

Anyway, once the script goes from an angle of 180 to -170, the optimisation fails, returning “Maximum dynamic_level reached”.

If I grep the .dat file for (D(2,4,5,10)):

	[D(2,4,5,10)[180.00,180.00]     179.68495        0.00000        0.21580     179.90075
	[D(2,4,5,10)[180.00,180.00]=         3.139860         179.900748
	[D(2,4,5,10)[180.00,180.00]     179.90075        0.00000        0.06798     179.96873
	[D(2,4,5,10)[180.00,180.00]=         3.141047         179.968731
	[D(2,4,5,10)[-170.00,-170.00]=         3.139860         179.900748
	[D(2,4,5,10)[-170.00,-170.00]     179.90075        0.00000       75.11985     255.02060

It appears as if there is some error parsing the angle value, jumping the angle up to 255?

When I tried by iterating using only positive values (i.e. np.arange(0,360,10)), the same error occurred at an angle of 200, so 10 more than the previous attempt.

When I tried doing the iterations in reverse, the script ran as far as angle=40, after which the system never converged, returning:
“invalid value encountered in scalar power”.

Am I doing something horribly wrong?

Many thanks in advance,
Jackson

          Psi4: An Open-Source Ab Initio Electronic Structure Package
                               Psi4 1.9.1 release

                         Git: Rev {} zzzzzzz
2

This is an interesting calculation. I have reproduced your results, and I’m looking into it. The first hurdle was obviously related to the discontinuity in the dihedral value, but, as you discovered, going around the other way is also failing for a different reason.

3

Thanks a bunch. With respect to the dihedral discontinuity, what’s the best way to do a full 360 degree scan with 10 degree iterations? Since once I get to 200 degrees, I seem to get the maximum dynamic_level reached error as well…
Thanks again!!

4

The geomeTRIC optimizer may be able to complete your dihedral scan; I was able to complete the job locally. Note that you need to change the atom indices to zero-based in the input file (attached below).

Thanks,

  • Lee-Ping
molecule = psi4.geometry("""
24
CHA.3.xyz
 C   -3.443370463455    0.647769656644    0.523542950280
 C   -1.959792017625    0.803964008843    0.703160331616
 O   -1.440827881011    1.679894908141    1.314167916183
 O   -1.296886858140   -0.168443235859    0.090654887976
 C    0.132792874268   -0.186700385789    0.167190539733
 C    0.737221850414    0.717431911434   -0.900072828218
 C    2.263890848138    0.601938019365   -0.910288899292
 C    2.718587007872   -0.849309082851   -1.078870699323
 C    2.089559565477   -1.755022848729   -0.017683497460
 C    0.563692287603   -1.635031555035   -0.008107529092
 H   -3.960379432774    1.451157007695    1.040215076678
 H   -3.757967278367   -0.319461737706    0.916057660454
 H   -3.683766745564    0.666052536222   -0.539809313296
 H    0.422313924444    0.172445636592    1.153883993526
 H    0.432840937126    1.749073688730   -0.720655923361
 H    0.334130586384    0.419982319792   -1.873857459669
 H    2.662872021991    1.000338867476    0.029278427713
 H    2.675192114576    1.224756843754   -1.708071556679
 H    3.808639452287   -0.911760818497   -1.028868558274
 H    2.430318344225   -1.205356867369   -2.074522820986
 H    2.481343810016   -1.481953560380    0.968748324716
 H    2.376730678543   -2.795566543583   -0.186401470489
 H    0.132170590827   -2.241898250449    0.790828821700
 H    0.153168828402   -2.010315508120   -0.950552627170

""")


set {
    basis cc-pvdz
    scf_type df
    intrafrag_step_limit 0.25
    g_convergence gau
}

def define_angles(psi=0):
    angles = np.arange(0, 360, 10) - 170
    idx = np.where(angles == int(psi))[0][0]
    return(np.concatenate([angles[idx:], angles[:idx]]))

# C-dihedral has initial angle of 153, so we start with 150...
Cdih_range = define_angles(150)
PES_C = []
for psi in Cdih_range:

    geometric_keywords = {
      'coordsys' : 'tric',
      'constraints' : {
      'set' : [{'type'    : 'dihedral',
                'indices' : [1, 3, 4, 9],
                'value'   : psi }]
      }
    }
    E = optimize('scf', engine='geometric', optimizer_keywords=geometric_keywords)
    out = f"CHA.C.{psi}.xyz"
    print(psi, E)
    PES_C.append((psi, E))
    molecule.save_xyz_file(out,1)
5

Hi, it looks like this worked locally for me too! Thank you very much for the help. I will still have some work and learning to do in order to intrepret the results, but for now I managed to make it work.
Thank you very much!
Jackson

6

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.