import mlqm import numpy as np import psi4 au2debye = psi4.constants.dipmom_au2debye # make directory lists name = 'CH3OH' # name of data dir # build directories for PESs #savedir = 'diatomics/' + name + '/' # where to save files #bdir = savedir + name.lower() + '_pes' # relative path to data dir #b = 0.5 #pes = [b] #for i in range(0,199): # b += (2.0 - 0.5) / 200 # pes.append(b) #pes = np.round(pes,4) #dirs = [(bdir + '/{}'.format(dis)) for dis in pes] # build directories for small molecules savedir = 'solvents/' + name + '/' # where to save files #savedir = 'chiral/' + name + '/' # where to save files bdir = savedir + name.lower() # relative path to data dir vals = ["%03d" % i for i in range(1,151)] # directory numbers dirs = [bdir + '/' + val for val in vals] # total directories # setup for MBTRs n_pts = 150 # points for discretizing cut_type = "full" # type of cutoff for DTR cut_val = 1000 # value of cutoff for DTR (if cut_type != 'full') # grab all relevant data w/ grabber function print("Grabbing scalars and densities . . .") results = mlqm.datahelper.grabber( dirs, varnames=[ 'CCSD CORRELATION ENERGY', 'MP2 CORRELATION ENERGY', 'HF TOTAL ENERGY', 'SCF DIPOLE X', 'SCF DIPOLE Y', 'SCF DIPOLE Z', 'MP2 DIPOLE X', 'MP2 DIPOLE Y', 'MP2 DIPOLE Z', 'CC DIPOLE X', 'CC DIPOLE Y', 'CC DIPOLE Z', ], fnames=['mp2_opdm.npy']) # additional steps in dirs: # # harvest nuclear dipoles # # compute total, nuclear, and electronic dipole moments # # transform densities # # harvest amps from output.dat print("Harvesting additional output data . . .") dip = {} # total dipole (nuc + scf + cc) e_dip = {} # electronic dipole (scf + cc) mp2_dip = {} # mp2 dip (no scf, no nuc) cc_dip = {} # cc dip (no scf, no nuc) ao_opdm = {} amps = {} for d in dirs: # total dipole (Debye) tx = results['CC DIPOLE X'][d] ty = results['CC DIPOLE Y'][d] tz = results['CC DIPOLE Z'][d] dip[d] = (tx**2 + ty**2 + tz**2)**0.5 # nuclear dipole (given in au, converted to debye) ndip = psi4.core.Wavefunction.from_file(d+'/ccsd_wfn.npy').molecule().nuclear_dipole() nx = ndip[0] * au2debye ny = ndip[1] * au2debye nz = ndip[2] * au2debye # electronic dipole (debye) ex = tx - nx ey = ty - ny ez = tz - nz e_dip[d] = (ex**2 + ey**2 + ez**2)**0.5 # total SCF dipole (given in debye) sx = results['SCF DIPOLE X'][d] sy = results['SCF DIPOLE Y'][d] sz = results['SCF DIPOLE Z'][d] # CC dipole (given in debye) cx = tx - sx cy = ty - sy cz = tz - sz cc_dip[d] = (cx**2 + cy**2 + cz**2)**0.5 # MP2 dipole (Debye) mx = results['MP2 DIPOLE X'][d] - sx my = results['MP2 DIPOLE Y'][d] - sy mz = results['MP2 DIPOLE Z'][d] - sz mp2_dip[d] = (mx**2 + my**2 + mz**2)**0.5 # # AO-basis OPDM # C = psi4.core.Wavefunction.from_file(d+'/mp2_wfn.npy').Ca().to_array() # ao_opdm[d] = C @ results['mp2_opdm.npy'][d] @ C.T # # # harvest amps # amps[d] = mlqm.datahelper.harvest_amps('mp2',outfile=d+'/output.dat')['t2'] # ## generate representations #print("Generating TATRs . . .") #tatrs = [mlqm.repgen.make_tatr('MP2',amps[d],x=150,st=0.05) for d in dirs] #print("Generating DTRs . . .") #dtrs = [mlqm.repgen.make_dtr(results['mp2_opdm.npy'][d],x=n_pts,cut_type=cut_type,cut_val=cut_val) for d in dirs] #print("Generating AO-basis DTRs . . .") #ao_dtrs = [mlqm.repgen.make_dtr(ao_opdm[d],x=n_pts,cut_type=cut_type,cut_val=cut_val) for d in dirs] # save to file print("Saving to file . . .") np.save('{}/{}_scf.npy'.format(savedir,name),list(results['HF TOTAL ENERGY'].values())) np.save('{}/{}_mp2.npy'.format(savedir,name),list(results['MP2 CORRELATION ENERGY'].values())) np.save('{}/{}_ccsd.npy'.format(savedir,name),list(results['CCSD CORRELATION ENERGY'].values())) np.save('{}/{}_dip.npy'.format(savedir,name),list(dip.values())) np.save('{}/{}_e_dip.npy'.format(savedir,name),list(e_dip.values())) np.save('{}/{}_mp2_dip.npy'.format(savedir,name),list(mp2_dip.values())) np.save('{}/{}_cc_dip.npy'.format(savedir,name),list(cc_dip.values())) #np.save('{}/{}_tatrs.npy'.format(savedir,name),tatrs) #np.save('{}/{}_dtrs.npy'.format(savedir,name),dtrs) #np.save('{}/{}_ao_dtrs.npy'.format(savedir,name),dtrs)