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