fragment MO guess for PSI4

fmo.in

# C1/RHF example of a fragment MO guess for the water dimer
# 

molecule wat {                                                                                                                                                                                                                                                                 
symmetry c1                                                                                                                                                                                                                                                                    
no_reorient
O      -1.769142     -0.076181     -0.000000
H      -2.065745      0.837492      0.000000
H      -0.809034      0.001317     -0.000000
--
O       1.283899      0.088882      0.000000
H       1.580006     -0.425760      0.756361
H       1.580006     -0.425760     -0.756361
}

set basis 3-21G
method='scf'
set guess read
set maxiter 30
set fail_on_maxiter false
set scf_type df

wat.fix_com(True)
wat.fix_orientation(True)

print "== monomer A =="
molA=wat.extract_subsets(1)
A_e_scf, Awfn = energy(method, molecule=molA, return_wfn=True)
print A_e_scf

A_occ = np.array(Awfn.Ca_subset("AO", "OCC").to_array()) 
A_virt = np.array(Awfn.Ca_subset("AO", "VIR").to_array())
Anocc=A_occ.shape[1]
Aoccpi=Awfn.doccpi()
Amo=A_occ.shape[0]
Anvirt=A_virt.shape[1]
Aeigen=np.array( Awfn.epsilon_a().to_array() )


print "== monomer B =="

molB=wat.extract_subsets(2)

B_e_scf, Bwfn = energy(method, molecule=molB, return_wfn=True)
print B_e_scf

B_occ = np.array(Bwfn.Ca_subset("AO", "OCC").to_array())
B_virt = np.array(Bwfn.Ca_subset("AO", "VIR").to_array())
Boccpi=Bwfn.doccpi()
Bnocc=B_occ.shape[1]
Bnvirt=B_virt.shape[1]
Bmo=B_occ.shape[0]
Beigen=np.array( Bwfn.epsilon_a().to_array() )

print "= promo wfn =="
clean()
ABnmo=Anocc+Bnocc+Anvirt+Bnvirt
AB=np.zeros((ABnmo,ABnmo))
ABeigen=np.zeros(ABnmo)
ABwfn = core.Wavefunction.build(wat, core.get_global_option('BASIS'))


# merge MOs and eigenvalues
# eigenvalue merge not needed
cnt=0
for i in range(Anocc):
   ABeigen[cnt]=Aeigen[i]
   for j in range(Amo):
      AB[j][cnt]=A_occ[j][i]
   cnt+=1

for i in range(Bnocc):
   ABeigen[cnt]=Beigen[i]
   for j in range(Bmo):
      AB[j+Amo][cnt]=B_occ[j][i]
   cnt+=1

# not needed
#for i in range(Anvirt):
#   ABeigen[cnt]=Aeigen[Anocc+i]
#   for j in range(Amo):
#      AB[j][cnt]=A_virt[j][i]
#   cnt+=1
#
#for i in range(Bnvirt):
#   ABeigen[cnt]=Beigen[Bnocc+i]
#   for j in range(Bmo):
#      AB[j+Amo][cnt]=B_virt[j][i]
#   cnt+=1


# occupied AB MOs
ABnocc=Anocc+Bnocc
AB_occ=AB[:,0:ABnocc]

# build new wfn object and set MOs as guess
# only MO_occ are saved
ABocc_Ca=core.Matrix.from_array(AB_occ)
AB_Ca=core.Matrix.from_array(AB)



# if you want a molden file without occupations
#new_wfn = scf_wavefunction_factory("scf", ABwfn, core.get_option('SCF', 'REFERENCE'))
#new_wfn.Ca().copy(ABocc_Ca)
#molden(new_wfn,filename='api.molden')


core.set_local_option('SCF', 'GUESS', 'READ')

# Write out product wfn orbitals for MO guess
scf_molecule = core.get_active_molecule()
read_orbitals = core.get_option('SCF', 'GUESS') is "READ"

name="wat"
fname = os.path.split(os.path.abspath(core.get_writer_file_prefix(name)))[1]
psi_scratch = core.IOManager.shared_object().get_default_path()
read_filename = os.path.join(psi_scratch, fname + ".180.npz")

write_filename = os.path.join(psi_scratch, fname + ".180.npz")
data = {}
data.update(ABocc_Ca.np_write(None, prefix="Ca_occ"))
data.update(ABocc_Ca.np_write(None, prefix="Cb_occ"))

data["reference"] = core.get_option('SCF', 'REFERENCE')
data["symmetry"] = wat.schoenflies_symbol()
data["BasisSet"] = ABwfn.basisset().name()
#print data

print write_filename
np.savez(write_filename, **data)
extras.register_numpy_file(write_filename)


print " == psi4 SCF  =="
activate(wat)
scf_wfn = scf_helper('scf', post_scf=False)
print core.get_variable('CURRENT ENERGY')