1#!/usr/bin/env python 2# Copyright 2014-2018 The PySCF Developers. All Rights Reserved. 3# 4# Licensed under the Apache License, Version 2.0 (the "License"); 5# you may not use this file except in compliance with the License. 6# You may obtain a copy of the License at 7# 8# http://www.apache.org/licenses/LICENSE-2.0 9# 10# Unless required by applicable law or agreed to in writing, software 11# distributed under the License is distributed on an "AS IS" BASIS, 12# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13# See the License for the specific language governing permissions and 14# limitations under the License. 15 16''' 17FCIDUMP functions (write, read) for real Hamiltonian 18''' 19 20import re 21from functools import reduce 22import numpy 23from pyscf import gto 24from pyscf import scf 25from pyscf import ao2mo 26from pyscf import symm 27from pyscf import __config__ 28 29DEFAULT_FLOAT_FORMAT = getattr(__config__, 'fcidump_float_format', ' %.16g') 30TOL = getattr(__config__, 'fcidump_write_tol', 1e-15) 31 32MOLPRO_ORBSYM = getattr(__config__, 'fcidump_molpro_orbsym', False) 33 34# Mapping Pyscf symmetry numbering to Molpro symmetry numbering for each irrep. 35# See also pyscf.symm.param.IRREP_ID_TABLE 36# https://www.molpro.net/info/current/doc/manual/node36.html 37ORBSYM_MAP = { 38 'D2h': (1, # Ag 39 4, # B1g 40 6, # B2g 41 7, # B3g 42 8, # Au 43 5, # B1u 44 3, # B2u 45 2), # B3u 46 'C2v': (1, # A1 47 4, # A2 48 2, # B1 49 3), # B2 50 'C2h': (1, # Ag 51 4, # Bg 52 2, # Au 53 3), # Bu 54 'D2' : (1, # A 55 4, # B1 56 3, # B2 57 2), # B3 58 'Cs' : (1, # A' 59 2), # A" 60 'C2' : (1, # A 61 2), # B 62 'Ci' : (1, # Ag 63 2), # Au 64 'C1' : (1,) 65} 66 67def write_head(fout, nmo, nelec, ms=0, orbsym=None): 68 if not isinstance(nelec, (int, numpy.number)): 69 ms = abs(nelec[0] - nelec[1]) 70 nelec = nelec[0] + nelec[1] 71 fout.write(' &FCI NORB=%4d,NELEC=%2d,MS2=%d,\n' % (nmo, nelec, ms)) 72 if orbsym is not None and len(orbsym) > 0: 73 fout.write(' ORBSYM=%s\n' % ','.join([str(x) for x in orbsym])) 74 else: 75 fout.write(' ORBSYM=%s\n' % ('1,' * nmo)) 76 fout.write(' ISYM=1,\n') 77 fout.write(' &END\n') 78 79 80def write_eri(fout, eri, nmo, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT): 81 npair = nmo*(nmo+1)//2 82 output_format = float_format + ' %4d %4d %4d %4d\n' 83 if eri.size == nmo**4: 84 eri = ao2mo.restore(8, eri, nmo) 85 86 if eri.ndim == 2: # 4-fold symmetry 87 assert(eri.size == npair**2) 88 ij = 0 89 for i in range(nmo): 90 for j in range(0, i+1): 91 kl = 0 92 for k in range(0, nmo): 93 for l in range(0, k+1): 94 if abs(eri[ij,kl]) > tol: 95 fout.write(output_format % (eri[ij,kl], i+1, j+1, k+1, l+1)) 96 kl += 1 97 ij += 1 98 else: # 8-fold symmetry 99 assert(eri.size == npair*(npair+1)//2) 100 ij = 0 101 ijkl = 0 102 for i in range(nmo): 103 for j in range(0, i+1): 104 kl = 0 105 for k in range(0, i+1): 106 for l in range(0, k+1): 107 if ij >= kl: 108 if abs(eri[ijkl]) > tol: 109 fout.write(output_format % (eri[ijkl], i+1, j+1, k+1, l+1)) 110 ijkl += 1 111 kl += 1 112 ij += 1 113 114def write_hcore(fout, h, nmo, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT): 115 h = h.reshape(nmo,nmo) 116 output_format = float_format + ' %4d %4d 0 0\n' 117 for i in range(nmo): 118 for j in range(0, i+1): 119 if abs(h[i,j]) > tol: 120 fout.write(output_format % (h[i,j], i+1, j+1)) 121 122 123def from_chkfile(filename, chkfile, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT, 124 molpro_orbsym=MOLPRO_ORBSYM, orbsym=None): 125 '''Read SCF results from PySCF chkfile and transform 1-electron, 126 2-electron integrals using the SCF orbitals. The transformed integrals is 127 written to FCIDUMP 128 129 Kwargs: 130 molpro_orbsym (bool): Whether to dump the orbsym in Molpro orbsym 131 convention as documented in 132 https://www.molpro.net/info/current/doc/manual/node36.html 133 ''' 134 mol, scf_rec = scf.chkfile.load_scf(chkfile) 135 mo_coeff = numpy.array(scf_rec['mo_coeff']) 136 nmo = mo_coeff.shape[1] 137 138 s = reduce(numpy.dot, (mo_coeff.conj().T, mol.intor_symmetric('int1e_ovlp'), mo_coeff)) 139 if abs(s - numpy.eye(nmo)).max() > 1e-6: 140 # Not support the chkfile from pbc calculation 141 raise RuntimeError('Non-orthogonal orbitals found in chkfile') 142 143 if mol.symmetry: 144 orbsym = symm.label_orb_symm(mol, mol.irrep_id, 145 mol.symm_orb, mo_coeff, check=False) 146 from_mo(mol, filename, mo_coeff, orbsym=orbsym, tol=tol, 147 float_format=float_format, molpro_orbsym=molpro_orbsym, 148 ms=mol.spin) 149 150def from_integrals(filename, h1e, h2e, nmo, nelec, nuc=0, ms=0, orbsym=None, 151 tol=TOL, float_format=DEFAULT_FLOAT_FORMAT): 152 '''Convert the given 1-electron and 2-electron integrals to FCIDUMP format''' 153 with open(filename, 'w') as fout: 154 write_head(fout, nmo, nelec, ms, orbsym) 155 write_eri(fout, h2e, nmo, tol=tol, float_format=float_format) 156 write_hcore(fout, h1e, nmo, tol=tol, float_format=float_format) 157 output_format = float_format + ' 0 0 0 0\n' 158 fout.write(output_format % nuc) 159 160def from_mo(mol, filename, mo_coeff, orbsym=None, 161 tol=TOL, float_format=DEFAULT_FLOAT_FORMAT, 162 molpro_orbsym=MOLPRO_ORBSYM, ms=0): 163 '''Use the given MOs to transfrom the 1-electron and 2-electron integrals 164 then dump them to FCIDUMP. 165 166 Kwargs: 167 molpro_orbsym (bool): Whether to dump the orbsym in Molpro orbsym 168 convention as documented in 169 https://www.molpro.net/info/current/doc/manual/node36.html 170 ''' 171 if getattr(mol, '_mesh', None): 172 raise NotImplementedError('PBC system') 173 174 if orbsym is None: 175 orbsym = getattr(mo_coeff, 'orbsym', None) 176 if molpro_orbsym and orbsym is not None: 177 orbsym = [ORBSYM_MAP[mol.groupname][i] for i in orbsym] 178 h1ao = scf.hf.get_hcore(mol) 179 h1e = reduce(numpy.dot, (mo_coeff.T, h1ao, mo_coeff)) 180 eri = ao2mo.full(mol, mo_coeff, verbose=0) 181 nuc = mol.energy_nuc() 182 from_integrals(filename, h1e, eri, h1e.shape[0], mol.nelec, nuc, ms, orbsym, 183 tol, float_format) 184 185def from_scf(mf, filename, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT, 186 molpro_orbsym=MOLPRO_ORBSYM): 187 '''Use the given SCF object to transfrom the 1-electron and 2-electron 188 integrals then dump them to FCIDUMP. 189 190 Kwargs: 191 molpro_orbsym (bool): Whether to dump the orbsym in Molpro orbsym 192 convention as documented in 193 https://www.molpro.net/info/current/doc/manual/node36.html 194 ''' 195 mol = mf.mol 196 mo_coeff = mf.mo_coeff 197 assert mo_coeff.dtype == numpy.double 198 199 h1e = reduce(numpy.dot, (mo_coeff.T, mf.get_hcore(), mo_coeff)) 200 if mf._eri is None: 201 if getattr(mf, 'exxdiv', None): # PBC system 202 eri = mf.with_df.ao2mo(mo_coeff) 203 else: 204 eri = ao2mo.full(mf.mol, mo_coeff) 205 else: # Handle cached integrals or customized systems 206 eri = ao2mo.full(mf._eri, mo_coeff) 207 orbsym = getattr(mo_coeff, 'orbsym', None) 208 if molpro_orbsym and orbsym is not None: 209 orbsym = [ORBSYM_MAP[mol.groupname][i] for i in orbsym] 210 nuc = mf.energy_nuc() 211 from_integrals(filename, h1e, eri, h1e.shape[0], mf.mol.nelec, nuc, 0, orbsym, 212 tol, float_format) 213 214 215def read(filename, molpro_orbsym=MOLPRO_ORBSYM): 216 '''Parse FCIDUMP. Return a dictionary to hold the integrals and 217 parameters with keys: H1, H2, ECORE, NORB, NELEC, MS, ORBSYM, ISYM 218 219 Kwargs: 220 molpro_orbsym (bool): Whether the orbsym in the FCIDUMP file is in 221 Molpro orbsym convention as documented in 222 https://www.molpro.net/info/current/doc/manual/node36.html 223 In return, orbsym is converted to pyscf symmetry convention 224 ''' 225 print('Parsing %s' % filename) 226 finp = open(filename, 'r') 227 228 data = [] 229 for i in range(10): 230 line = finp.readline().upper() 231 data.append(line) 232 if '&END' in line: 233 break 234 else: 235 raise RuntimeError('Problematic FCIDUMP header') 236 237 result = {} 238 tokens = ','.join(data).replace('&FCI', '').replace('&END', '') 239 tokens = tokens.replace(' ', '').replace('\n', '').replace(',,', ',') 240 for token in re.split(',(?=[a-zA-Z])', tokens): 241 key, val = token.split('=') 242 if key in ('NORB', 'NELEC', 'MS2', 'ISYM'): 243 result[key] = int(val.replace(',', '')) 244 elif key in ('ORBSYM',): 245 result[key] = [int(x) for x in val.replace(',', ' ').split()] 246 else: 247 result[key] = val 248 249 # Convert to molpr orbsym convert_orbsym 250 if 'ORBSYM' in result: 251 if molpro_orbsym: 252 # Guess which point group the orbsym belongs to. FCIDUMP does not 253 # save the point group information, the guess might be wrong if 254 # the high symmetry numbering of orbitals are not presented. 255 orbsym = result['ORBSYM'] 256 if max(orbsym) > 4: 257 result['ORBSYM'] = [ORBSYM_MAP['D2h'].index(i) for i in orbsym] 258 elif max(orbsym) > 2: 259 # Fortunately, without molecular orientation, B2 and B3 in D2 260 # are not distinguishable 261 result['ORBSYM'] = [ORBSYM_MAP['C2v'].index(i) for i in orbsym] 262 elif max(orbsym) == 2: 263 result['ORBSYM'] = [i-1 for i in orbsym] 264 elif max(orbsym) == 1: 265 result['ORBSYM'] = [0] * len(orbsym) 266 else: 267 raise RuntimeError('Unknown orbsym') 268 elif max(result['ORBSYM']) >= 8: 269 raise RuntimeError('Unknown orbsym convention') 270 271 norb = result['NORB'] 272 norb_pair = norb * (norb+1) // 2 273 h1e = numpy.zeros((norb,norb)) 274 h2e = numpy.zeros(norb_pair*(norb_pair+1)//2) 275 dat = finp.readline().split() 276 while dat: 277 i, j, k, l = [int(x) for x in dat[1:5]] 278 if k != 0: 279 if i >= j: 280 ij = i * (i-1) // 2 + j-1 281 else: 282 ij = j * (j-1) // 2 + i-1 283 if k >= l: 284 kl = k * (k-1) // 2 + l-1 285 else: 286 kl = l * (l-1) // 2 + k-1 287 if ij >= kl: 288 h2e[ij*(ij+1)//2+kl] = float(dat[0]) 289 else: 290 h2e[kl*(kl+1)//2+ij] = float(dat[0]) 291 elif k == 0: 292 if j != 0: 293 h1e[i-1,j-1] = float(dat[0]) 294 else: 295 result['ECORE'] = float(dat[0]) 296 dat = finp.readline().split() 297 298 idx, idy = numpy.tril_indices(norb, -1) 299 if numpy.linalg.norm(h1e[idy,idx]) == 0: 300 h1e[idy,idx] = h1e[idx,idy] 301 elif numpy.linalg.norm(h1e[idx,idy]) == 0: 302 h1e[idx,idy] = h1e[idy,idx] 303 result['H1'] = h1e 304 result['H2'] = h2e 305 finp.close() 306 return result 307 308def to_scf(filename, molpro_orbsym=MOLPRO_ORBSYM, mf=None, **kwargs): 309 '''Use the Hamiltonians defined by FCIDUMP to build an SCF object''' 310 ctx = read(filename, molpro_orbsym) 311 mol = gto.M() 312 mol.nelectron = ctx['NELEC'] 313 mol.spin = ctx['MS2'] 314 norb = mol.nao = ctx['NORB'] 315 if 'ECORE' in ctx: 316 mol.energy_nuc = lambda *args: ctx['ECORE'] 317 mol.incore_anyway = True 318 319 if 'ORBSYM' in ctx: 320 mol.symmetry = True 321 mol.groupname = 'N/A' 322 orbsym = numpy.asarray(ctx['ORBSYM']) 323 mol.irrep_id = list(set(orbsym)) 324 mol.irrep_name = [('IR%d' % ir) for ir in mol.irrep_id] 325 so = numpy.eye(norb) 326 mol.symm_orb = [] 327 for ir in mol.irrep_id: 328 mol.symm_orb.append(so[:,orbsym==ir]) 329 330 if mf is None: 331 mf = mol.RHF(**kwargs) 332 else: 333 mf.mol = mol 334 h1 = ctx['H1'] 335 idx, idy = numpy.tril_indices(norb, -1) 336 if h1[idx,idy].max() == 0: 337 h1[idx,idy] = h1[idy,idx] 338 else: 339 h1[idy,idx] = h1[idx,idy] 340 mf.get_hcore = lambda *args: h1 341 mf.get_ovlp = lambda *args: numpy.eye(norb) 342 mf._eri = ctx['H2'] 343 344 return mf 345 346def scf_from_fcidump(mf, filename, molpro_orbsym=MOLPRO_ORBSYM): 347 '''Update the SCF object with the quantities defined in FCIDUMP file''' 348 return to_scf(filename, molpro_orbsym, mf) 349 350scf.hf.SCF.from_fcidump = scf_from_fcidump 351 352if __name__ == '__main__': 353 import argparse 354 parser = argparse.ArgumentParser(description='Convert chkfile to FCIDUMP') 355 parser.add_argument('chkfile', help='pyscf chkfile') 356 parser.add_argument('fcidump', help='FCIDUMP file') 357 args = parser.parse_args() 358 359 # fcidump.py chkfile output 360 from_chkfile(args.fcidump, args.chkfile) 361