Hello, I am using psi4, version 1.9.1. I am trying to carry out a single point energy calculations, save the wfn, and then use that wfn to generate an ESP grid via the makeESPgrid function. I am including my script below:
def getEfieldBasis(mol_xyz, idx_atom1, idx_atom2):
#xyz, charge, spin, symmetry are arguments
mol = get_molecule(mol_xyz, 0, 1)
# Set the desired HF exchange percentage (e.g., 30% HF exchange)
hf_exchange_percentage = 0.20
#using the default params for b3lyp composition:
b = 0.72
c = 0.81
# Define the custom functional
custom_functional = {"name" : 'b3lyp_hx' + str(hf_exchange_percentage),
"x_functionals": {
"LDA_X": {"alpha": 0.10*(1 - hf_exchange_percentage)}, # 10% non-HF XC isLDA exchange
"GGA_X_B88": {"alpha": 0.90*(1 - hf_exchange_percentage)}, # 90% of non-HF XC is Becke 1988 exchange
},
"x_hf": {"alpha": hf_exchange_percentage }, # 20% Hartree-Fock exchange
"c_functionals": {
"LDA_C_VWN": {"alpha": 0.19}, # 19% VWN correlation
"GGA_C_LYP": {"alpha": 0.81}, # 81% Lee-Yang-Parr correlation
}
}
mybasis = 'def2-TZVP'
basis_lst = ['def2-sv(p)', 'def2-SVPD', 'def2-TZVP', 'def2-TZVPD']
for mybasis in basis_lst[0:4]:
psi4.set_options({'basis': mybasis, 'scf_type': 'df', "print": 2})
psi4.core.set_output_file('output.dat', False)
e, wfn = psi4.energy('scf', dft_functional=custom_functional, return_wfn=True)
wfn.to_file('./wfn.180.npy')
use_previous_wfn = True
#num_pts specifies the grid where 0 is a grid with onyl the points of interst themselves and out-lying points at
#same interval are specified
avg_Efield_lst = []
num_pts_lst = []
times_lst = []
for num_pts in [4]:
print(f'Number of pts: {num_pts}')
start_time = time.time()
[Vvals, Exvals, Eyvals, Ezvals, average_Efield] = makeESPgrid(mol, num_pts, idx_atom1, idx_atom2, mybasis, use_previous_wfn)
num_pts_lst.append(num_pts)
avg_Efield_lst.append(average_Efield)
end_time = time.time()
times_lst.append(end_time - start_time)
df = pd.DataFrame({'Num Grid Pts' : num_pts_lst, 'avgEfield': avg_Efield_lst, 'Time': times_lst })
csv_name = 'HFx' + str(hf_exchange_percentage) + '_'+ mybasis + '_gridptsEfield.csv'
df.to_csv(csv_name)
idx_atom1 = 0
idx_atom2 = 1
mol_xyz = 'input.xyz'
getEfieldBasis(mol_xyz, idx_atom1, idx_atom2)
I am running these calculations for an acetone molecule in a solvated environment (with varying sizes of solvation shell). These calculations work without a problem when running below about 65 atoms, but crash before even finishing the B3LYP/def2-SVP single point above that number of atoms. The crash always happens right after “==> Integral Setup <==” is written and I can find no warnings or error messages. Note that in my jobscript, I explicitly change the scratch directory “export PSI_SCRATCH=‘./’” to avoid running into scratch space issues and I allocated 40GB to the job which should be more than enough for this calculation. I have checked local memory and scratch memory and both should be sufficient.