Commit f336ad1b authored by hebrewsnabla's avatar hebrewsnabla
Browse files

add frag_sol

parent c0476c8a
import frag_sol
xyz = '''H 0.0 0.0 0.0
O 0.0 0.0 1.0
O 0.0 0.0 2.0
H 0.0 0.0 3.0
'''
#bas = 'cc-pvdz'
gjfhead = '''%nproc=1
%mem=1gb
#p b3lyp/cc-pvdz scrf=read out=wfn
title
'''
scrfhead = '''eps=1.0
epsinf=1.0
'''
wfnpath = 'D:\\xxx\\yyy\\frag'
## the optput wfn will be frag0.wfn, frag1.wfn, ...
frag_sol.from_frag(xyz, [[0,1],[2,3]], [0,0],[1,-1], gjfhead, scrfhead, 'frag', wfnpath)
from pyscf import gto
import radii
def from_frag(xyz, frags, chgs, spins, gjfhead='', scrfhead='', gjfname='', wfnpath=''):
mol = gto.Mole()
mol.atom = xyz
# mol.basis = bas
mol.verbose = 1
mol.build()
guess_frag(mol, frags, chgs, spins, gjfhead, scrfhead, gjfname, wfnpath)
def spin_p2g(spin):
if spin>0:
spin = spin + 1
elif spin < 0:
spin = spin - 1
return spin
def guess_frag(mol, frags, chgs, spins, gjfhead, scrfhead, gjfname, wfnpath):
'''
frags: e.g. [[0], [1]] for N2
chgs: e.g. [0, 0] for N2
spins: e.g. [3, -3] for N2
'''
#mol.build()
print('**** generating fragments ****')
atom = mol.format_atom(mol.atom, unit=1)
#print(atom)
#fraga, fragb = frags
#chga, chgb = chgs
#spina, spinb = spins
allatom = range(len(atom))
for k in range(len(frags)):
frag = frags[k]
chg = chgs[k]
spin = spins[k]
g_spin = spin_p2g(spin)
atomk = [atom[i] for i in frag]
atomother = [atom[i] for i in allatom if i not in frag]
print('fragment %d, chg %d, spin %d' % (k, chg, spin))
#print(atomk)
with open(gjfname+'%d.gjf'%k, 'w') as f:
f.write(gjfhead)
f.write('%d %d\n' % (chg, g_spin))
for a in atomk:
f.write('%s %10.5f %10.5f %10.5f\n' % (a[0], a[1][0], a[1][1], a[1][2]))
f.write('\n')
f.write(scrfhead)
f.write('ExtraSph=%d\n\n' % len(atomother))
for b in atomother:
rad = radii.uff_radii[b[0]] / 2.0
f.write(' %10.5f %10.5f %10.5f %10.5f\n' % (b[1][0], b[1][1], b[1][2], rad))
f.write('\n')
f.write(wfnpath + '%d.wfn'%k + '\n')
f.write('\n')
#import hbonda
def smd_oxygen_rad(solvent):
if solvent in hbonda.hbonda:
alpha = hbonda.hbonda[solvent]
else:
alpha = 0.0
if alpha>0.429:
return 1.52
else:
return 1.52 + 1.8*(0.43-alpha)
smd_radii = {
# fit by Trular
'H': 1.2,
'C': 1.85,
'N': 1.89,
#'O': '--',
'F' : 1.73,
'Si': 2.47,
'P' : 2.12,
'S' : 2.49,
'Cl': 2.38,
'Br': 3.06,
# other: van der Waals radius of Bondi
'He': 1.40,
'Li': 1.82,
'Be': 1.45,
'B ': 1.80,
'Ne': 1.54,
'Na': 2.27,
'Mg': 1.73,
'Al': 2.30,
'Ar': 1.88,
'K' : 2.75,
#'Ca': '--',
'Ga': 1.87,
'Ge': 2.19,
'As': 1.85,
'Se': 1.90,
'Kr': 2.02,
'In': 1.93,
'Sn': 2.17,
#'Sb': '--',
'Te': 2.06,
'I' : 1.98,
'Xe': 2.16,
'Tl': 1.96,
'Pb': 2.02,
'Ni': 1.63,
'Cu': 1.40,
'Zn': 1.39,
'Pd': 1.63,
'Ag': 1.72,
'Cd': 1.62,
'Pt': 1.75, # 1.70-1.80?
'Au': 1.66,
'Hg': 1.70,
'U' : 1.86
}
uff_radii = {
# 2*rad
'H' : 2.886e+00,
'He': 2.362e+00,
'Li': 2.451e+00,
'Be': 2.745e+00,
'B' : 4.083e+00,
'C' : 3.851e+00,
'N' : 3.660e+00,
'O' : 3.500e+00,
'F' : 3.364e+00,
'Ne': 3.243e+00,
'Na': 2.983e+00,
'Mg': 3.021e+00,
'Al': 4.499e+00,
'Si': 4.295e+00,
'P' : 4.147e+00,
'S' : 4.035e+00,
'Cl': 3.947e+00,
'Ar': 3.868e+00,
'K' : 3.812e+00,
'Ca': 3.399e+00,
'Sc': 3.295e+00,
'Ti': 3.175e+00,
'V' : 3.144e+00,
'Cr': 3.023e+00,
'Mn': 2.961e+00,
'Fe': 2.912e+00,
'Co': 2.872e+00,
'Ni': 2.834e+00,
'Cu': 3.495e+00,
'Zn': 2.763e+00,
'Ga': 4.383e+00,
'Ge': 4.280e+00,
'As': 4.230e+00,
'Se': 4.205e+00,
'Br': 4.189e+00,
'Kr': 4.141e+00,
'Rb': 4.114e+00,
'Sr': #(+2)
3.641e+00,
'Y' : #(+3)
3.345e+00,
'Zr': #(+4)
3.124e+00,
'Nb': #(+5)
3.165e+00,
'Mo': #(+6)
3.052e+00,
'Tc': #(+5)
2.998e+00,
'Ru': #(+2)
2.963e+00,
'Rh': #(+3)
2.929e+00,
'Pd': #(+2)
2.899e+00,
'Ag': #(+1)
3.148e+00,
'Cd': #(+2)
2.848e+00,
'In': 4.463e+00,
'Sn': 4.392e+00,
'Sb': 4.420e+00,
'Te': 4.470e+00,
'I' : 4.50e+00,
'Xe': 4.404e+00,
'Cs': 4.517e+00,
'Ba': #(+2)
3.703e+00,
'La': #(+3)
3.522e+00,
'Ce': #(+3)
3.556e+00,
'Pr': #(+3)
3.606e+00,
'Nd': #(+3)
3.575e+00,
'Pm': #(+3)
3.547e+00,
'Sm': # (+3)
3.520e+00,
'Eu': # (+3)
3.493e+00,
'Gd': # (+3)
3.368e+00,
'Tb': # (+3)
3.451e+00,
'Dy': # (+3)
3.428e+00,
'Ho': # (+3)
3.409e+00,
'Er': # (+3)
3.391e+00,
'Tm': # (+3)
3.374e+00,
'Yb': # (+3)
3.355e+00,
'Lu': # (+3)
3.640e+00,
'Hf': # (+4)
3.141e+00,
'Ta': # (+5)
3.170e+00,
'W' : #(+4,+6)
3.069e+00,
'Re': # (+5,+7)
2.954e+00,
'Os': # (+6)
3.120e+00,
'Ir': # (+3)
2.840e+00,
'Pt': 2.754e+00,
'Au': 3.293e+00,
'Hg': 2.705e+00,
'Tl': 4.347e+00,
'Pb': 4.297e+00,
'Bi': # (+3)
4.370e+00,
'Po': # (+2)
4.709e+00,
'At': 4.750e+00,
'Rn': # (+4)
4.765e+00,
'Fr':
4.900e+00,
'Ra': # (+2)
3.677e+00,
'Ac': # (+3)
3.478e+00,
'Th': # (+4)
3.396e+00,
'Pa': # (+4)
3.424e+00,
'U': #(+4)
3.395e+00,
'Np': # (+4)
3.424e+00,
'Pu': # (+4)
3.424e+00,
'Am': # (+4)
3.381e+00,
'Cm': # (+3)
3.326e+00,
'Bk': # (+3)
3.339e+00,
'Cf': # (+3)
3.313e+00,
'Es': # (+3)
3.299e+00,
'Fm': # (+3)
3.286e+00,
'Md': # (+3)
3.274e+00,
'No': # (+3)
3.248e+00,
'Lw': #(+3)
3.236e+00,
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment