Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# bar{C}_{pi}^1 ~= C_{pi}^1 -
if sscobj.mb.upper() == 'RMB':
orbo = mo_coeff[:,mo_occ> 0]
orbv = mo_coeff[:,mo_occ==0]
n4c = mo_coeff.shape[0]
n2c = n4c // 2
c = lib.param.LIGHT_SPEED
orbvS_T = orbv[n2c:].conj().T
for ia in atmlst:
mol.set_rinv_origin(mol.atom_coord(ia))
a01int = mol.intor('int1e_sa01sp_spinor', 3)
for k in range(3):
s1 = orbvS_T.dot(a01int[k].conj().T).dot(orbo[n2c:])
mo1[ia*3+k] -= s1 * (.25/c**2)
logger.timer(sscobj, 'solving mo1 eqn', *cput1)
return mo1, mo_e1
def get_pnucp(mydf, kpts=None):
cell = mydf.cell
if kpts is None:
kpts_lst = numpy.zeros((1,3))
else:
kpts_lst = numpy.reshape(kpts, (-1,3))
log = logger.Logger(mydf.stdout, mydf.verbose)
t1 = t0 = (time.clock(), time.time())
nkpts = len(kpts_lst)
nao = cell.nao_nr()
nao_pair = nao * (nao+1) // 2
Gv, Gvbase, kws = cell.get_Gv_weights(mydf.mesh)
charge = -cell.atom_charges()
kpt_allow = numpy.zeros(3)
coulG = tools.get_coulG(cell, kpt_allow, mesh=mydf.mesh, Gv=Gv)
coulG *= kws
if mydf.eta == 0:
wj = numpy.zeros((nkpts,nao_pair), dtype=numpy.complex128)
SI = cell.get_SI(Gv)
vG = numpy.einsum('i,ix->x', charge, SI) * coulG
if cell.dimension == 1 or cell.dimension == 2:
_basis = mol._basis
auxbasis = {}
for k in _basis:
if isinstance(_basis[k], str):
balias = gto.basis._format_basis_name(_basis[k])
if gto.basis._is_pople_basis(balias):
balias = balias.split('g')[0] + 'g'
if balias in DEFAULT_AUXBASIS:
if mp2fit:
auxb = DEFAULT_AUXBASIS[balias][1]
else:
auxb = DEFAULT_AUXBASIS[balias][0]
if auxb is not None and gto.basis.load(auxb, k):
auxbasis[k] = auxb
logger.info(mol, 'Default auxbasis %s is used for %s %s',
auxb, k, _basis[k])
if len(auxbasis) != len(_basis):
# Some AO basis not found in DEFAULT_AUXBASIS
auxbasis, auxdefault = aug_etb(mol), auxbasis
auxbasis.update(auxdefault)
aux_etb = set(auxbasis) - set(auxdefault)
if aux_etb:
logger.info(mol, 'Even tempered Gaussians are generated as '
'DF auxbasis for %s', ' '.join(aux_etb))
for k in aux_etb:
logger.debug(mol, ' ETB auxbasis for %s %s', k, auxbasis[k])
return auxbasis
def dump_flags(self, verbose=None):
log = logger.new_logger(self, verbose)
log.info('\n')
log.info('******** %s for %s ********',
self.__class__, self._scf.__class__)
if self.singlet:
log.info('nstates = %d singlet', self.nstates)
else:
log.info('nstates = %d triplet', self.nstates)
log.info('wfnsym = %s', self.wfnsym)
log.info('conv_tol = %g', self.conv_tol)
log.info('eigh lindep = %g', self.lindep)
log.info('eigh level_shift = %g', self.level_shift)
log.info('eigh max_space = %d', self.max_space)
log.info('eigh max_cycle = %d', self.max_cycle)
log.info('chkfile = %s', self.chkfile)
log.info('max_memory %d MB (current use %d MB)',
self.max_memory, lib.current_memory()[0])
def run_mfs(mf, mols_a, mols_b):
'''perform a set of calculations on given two sets of molecules'''
nconfigs = len(mols_a)
dm0 = mf.make_rdm1()
mflist = []
for i in range(nconfigs):
mf1 = copy_mf(mf, mols_a[i])
mf2 = copy_mf(mf, mols_b[i])
mf1.kernel(dm0=dm0)
mf2.kernel(dm0=dm0)
if not (mf1.converged):
logger.warn(mf, "%ith config mf1 not converged", i)
if not (mf2.converged):
logger.warn(mf, "%ith config mf2 not converged", i)
mflist.append((mf1, mf2))
return mflist
def ecp_int(cell, kpts=None):
from pyscf.pbc.df import incore
if kpts is None:
kpts_lst = numpy.zeros((1,3))
else:
kpts_lst = numpy.reshape(kpts, (-1,3))
cell, contr_coeff = gto.cell._split_basis(cell)
lib.logger.debug1(cell, 'nao %d -> nao %d', *(contr_coeff.shape))
ecpcell = gto.Cell()
ecpcell._atm = cell._atm
# append a fictitious s function to mimic the auxiliary index in pbc.incore.
# ptr2last_env_idx to force PBCnr3c_fill_* function to copy the entire "env"
ptr2last_env_idx = len(cell._env) - 1
ecpbas = numpy.vstack([[0, 0, 1, 1, 0, ptr2last_env_idx, 0, 0],
cell._ecpbas]).astype(numpy.int32)
ecpcell._bas = ecpbas
ecpcell._env = cell._env
# In pbc.incore _ecpbas is appended to two sets of cell._bas and the
# fictitious s function.
cell._env[AS_ECPBAS_OFFSET] = cell.nbas * 2 + 1
cell._env[AS_NECPBAS] = len(cell._ecpbas)
# shls_slice of auxiliary index (0,1) corresponds to the fictitious s function
shls_slice = (0, cell.nbas, 0, cell.nbas, 0, 1)
self.t2_L2iL = None
self.t2_L_group = None
self.t2_sym_group = None
'''
self.t3_iLs_idx = None
self.t3_iL2L = None
self.t3_L2iL = None
self.t3_L_group = None
self.t3_sym_group = None
'''
self.shltrip_cen_idx = get_shl_center_idx(cell, 3)
self.nL = len(Ls)
self.nL2 = self.nL * self.nL
lib.logger.debug(self, 'nL2 = %s', self.nL2)
nshltrip = len(self.shltrip_cen_idx)
self.buf = np.memmap('petite_list', dtype=np.int32, shape=(nshltrip,self.nL2,2), mode='w+')
if k == 0:
self.stdout.write('%-15.12g ' % x[0])
else:
self.stdout.write(' '*32+'%-15.12g ' % x[0])
for c in x[1:]:
self.stdout.write(' %4.12g' % c)
self.stdout.write('\n')
if self.verbose >= logger.INFO:
self.stdout.write('\n')
logger.info(self, 'nuclear repulsion = %.15g', self.energy_nuc())
if self.symmetry:
if self.topgroup == self.groupname:
logger.info(self, 'point group symmetry = %s', self.topgroup)
else:
logger.info(self, 'point group symmetry = %s, use subgroup %s',
self.topgroup, self.groupname)
for ir in range(self.symm_orb.__len__()):
logger.info(self, 'num. orbitals of irrep %s = %d',
self.irrep_name[ir], self.symm_orb[ir].shape[1])
logger.info(self, 'number of shells = %d', self.nbas)
logger.info(self, 'number of NR pGTOs = %d', self.npgto_nr())
logger.info(self, 'number of NR cGTOs = %d', self.nao_nr())
logger.info(self, 'basis = %s', self.basis)
logger.info(self, 'ecp = %s', self.ecp)
if self.verbose >= logger.DEBUG2:
for i in range(len(self._bas)):
exps = self.bas_exp(i)
logger.debug1(self, 'bas %d, expnt(s) = %s', i, str(exps))
logger.info(self, 'CPU time: %12.2f', time.clock())
return self
def grad(self, mo_energy=None, mo_coeff=None, mo_occ=None, atmlst=None):
cput0 = (time.clock(), time.time())
if mo_energy is None: mo_energy = self._scf.mo_energy
if mo_coeff is None: mo_coeff = self._scf.mo_coeff
if mo_occ is None: mo_occ = self._scf.mo_occ
if atmlst is None:
atmlst = range(self.mol.natm)
if self.verbose >= logger.WARN:
self.check_sanity()
if self.verbose >= logger.INFO:
self.dump_flags()
de = self.grad_elec(mo_energy, mo_coeff, mo_occ, atmlst)
self.de = de = de + self.grad_nuc(atmlst=atmlst)
logger.note(self, '--------------')
logger.note(self, ' x y z')
for k, ia in enumerate(atmlst):
logger.note(self, '%d %s %15.9f %15.9f %15.9f', ia,
self.mol.atom_symbol(ia), de[k,0], de[k,1], de[k,2])
logger.note(self, '--------------')
logger.timer(self, 'SCF gradients', *cput0)
return self.de
mo_ea = mo_energy.mo_ea
mo_eb = mo_energy.mo_eb
else:
mo_ea = mo_eb = mo_energy
nmo = mo_ea.size
mo_occ = numpy.zeros(nmo)
if getattr(mf, 'nelec', None) is None:
nelec = mf.mol.nelec
else:
nelec = mf.nelec
ncore = nelec[1]
nocc = nelec[0]
nopen = abs(nocc - ncore)
mo_occ = _fill_rohf_occ(mo_energy, mo_ea, mo_eb, ncore, nopen)
if mf.verbose >= logger.INFO and nocc < nmo and ncore > 0:
ehomo = max(mo_energy[mo_occ> 0])
elumo = min(mo_energy[mo_occ==0])
if ehomo+1e-3 > elumo:
logger.warn(mf, 'HOMO %.15g >= LUMO %.15g', ehomo, elumo)
else:
logger.info(mf, ' HOMO = %.15g LUMO = %.15g', ehomo, elumo)
if nopen > 0 and mf.verbose >= logger.DEBUG:
core_idx = mo_occ == 2
open_idx = mo_occ == 1
vir_idx = mo_occ == 0
logger.debug(mf, ' Roothaan | alpha | beta')
logger.debug(mf, ' Highest 2-occ = %18.15g | %18.15g | %18.15g',
max(mo_energy[core_idx]),
max(mo_ea[core_idx]), max(mo_eb[core_idx]))
logger.debug(mf, ' Lowest 0-occ = %18.15g | %18.15g | %18.15g',
min(mo_energy[vir_idx]),