4
$\begingroup$

I was using psi4 to run vibrational frequency analysis with small molecules (in a jupyter notebook) and it works just fine with molecules like methanol and water but once I move up to ethanol or acetone, I can see the prompt can do the computations just fine but it gets stuck on saving the results file to my notebook (the prompt just says 'Saving file at x.ipynb'). Has anyone run into this error before? If so, how did you fix it?

code:

#import psi4 and set memory

import psi4

psi4.set_memory('10 GB')

from rdkit import Chem # main tools
from rdkit.Chem import AllChem # additional tools, including 3D

smiles = 'CCO' # smiles string input
mol = Chem.MolFromSmiles(smiles) # initialize molecule

mol_h = Chem.AddHs(mol) # adding explicit Hs for 3D generation
cid = AllChem.EmbedMolecule(mol_h) # returns the id of the generated conformer,
                                 # and -1 if no conformers were generated

AllChem.MMFFOptimizeMolecule(mol_h) # optimize molecule with MMFF94 forcefield
mol_xyz=Chem.rdmolfiles.MolToXYZBlock(mol_h) # save to xyz string

print(mol_xyz)

#check structure of molecule
mol

#convert xyz file to psi4 geometry file
qmol = psi4.driver.qcdb.Molecule.from_string(mol_xyz, dtype='xyz+')
mol_psi4 = psi4.geometry(qmol.create_psi4_string_from_molecule())

#initialize molecule 
mol_psi4.update_geometry()

#optimize geometry before calculating frequencies
psi4.set_options({'reference': 'rhf'})
psi4.optimize('scf/cc-pvdz', dft_functional = 'b3lyp', molecule=mol_psi4)

#calculate vibrational frequencies
scf_e, scf_wfn = psi4.frequency('scf/cc-pvdz', molecule=mol_psi4, return_wfn=True)
print('done')

#save and print frequencies + intesity
freq=scf_wfn.frequency_analysis['omega'][2]
for i in freq:
    print(i)

print('')
print('Done!!!')
$\endgroup$
2
  • $\begingroup$ +1. When you say "it gets stuck on saving the results file to my notebook", what is being printed? Could you please show us an example of what gets printed when you do a calculation on H2O for example? $\endgroup$ Commented Apr 23, 2023 at 21:37
  • $\begingroup$ @meheret Does my answer solve your question? If yes, could you please mark the answer as accepted by clicking on the tick mark next to it $\endgroup$ Commented May 18, 2023 at 18:26

1 Answer 1

4
$\begingroup$

I ran the calculation for ethanol and it took me 0.78 mins.

System Specs:

OS: MacOS Monterey 12.5.1
2,6 GHz 6-Core Intel Core i7
32 GB 2667 MHz DDR4
IDE: PyCharm 2022.3.2 (Community Edition)

Software Specs:

Virtual Env: Conda env 
Conda version: 23.1.0
Python version: 3.10.10
PSI4 version:  1.6.1+5b9f6e3   
PSI4 build : py310h581b2b9_0

In the following code I am directly using the xyz string obtained from rdkit. (The reason is that it is a hassle to install both rdkit and psi4 on the same conda env on OSX. rdkit needs python version <=3.7 and psi4 seems to work well only on 3.10 onwards)

Code:

import time
s = time.time()
mol_psi4 = psi4.geometry('''9

C     -0.887089    0.175064   -0.012535
C      0.460489   -0.515516   -0.046535
O      1.442965    0.307267    0.565572
H     -0.847478    1.127768   -0.550817
H     -1.658782   -0.453327   -0.465842
H     -1.176944    0.403676    1.018306
H      0.768712   -0.724328   -1.075460
H      0.419486   -1.462073    0.500177
H      1.478640    1.141468    0.067135
''')

#optimize geometry before calculating frequencies
psi4.set_options({'reference': 'rhf'})
psi4.optimize('scf/cc-pvdz', dft_functional = 'b3lyp', molecule=mol_psi4)

#calculate vibrational frequencies
scf_e, scf_wfn = psi4.frequency('scf/cc-pvdz', molecule=mol_psi4, return_wfn=True)
print('done')

#save and print frequencies + intesity
freq=(scf_wfn.frequency_analysis['omega'].data)
for i in freq:
    print(i)
e = time.time()

print(f'time taken mins: {(e-s)/60}')
print('')
print('Done!!!')

Also here are the output frequencies:

7.507156968262905e-06j
(4.908897302913366e-06+0j)
(1.5197889898678715e-05+0j)
(2.9472406989666606e-05+0j)
(3.222627049809006e-05+0j)
(4.2094538779691455e-05+0j)
(401.17046333897605+0j)
(462.0977321879201+0j)
(688.388935344985+0j)
(908.0756382425054+0j)
(951.6293253661514+0j)
(1133.4970378368564+0j)
(1195.0895426135098+0j)
(1238.377261101374+0j)
(1403.076854517634+0j)
(1506.414762164892+0j)
(1537.260320704411+0j)
(1574.0259056577736+0j)
(1596.5805492658997+0j)
(1604.014633238946+0j)
(1637.0354072816572+0j)
(3036.6416405029995+0j)
(3081.9122987732476+0j)
(3124.28961913401+0j)
(3149.4983929335085+0j)
(3168.1054762068184+0j)
(3824.5089839751754+0j)
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .