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)