Single point job stops over 65 atoms; no error message

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.

I got the same thing in a 24-atom cluster because my geometry was passing through itself, causing a basis set construction error. I don’t think I can help you beyond that but check that your input is right and that you didn’t make any periodicity mistakes.