pyscf - Githubissues

Apologies for the delay, below I have pasted my scripts to calculate the La-H binding curve. It is a bit heavy-handed but I have tried to keep it as close to how I originally ran it. Hopefully this example is sufficient to show how I got the data. I first calculated the energy as a function of separation with minimal manipulation, then transformed it into tables that were used in the study.

pyscf script:

from pyscf import gto, scf, cc
import numpy as np

mol = gto.Mole()

mol.basis = {'La': gto.basis.load('path/to/def2-qzvppd.1.nw', 'La'),
        'H': 'cc-pv5z'}
mol.ecp = {'La': gto.basis.load_ecp('path/to/def2-qzvppd.1.nw', 'La')}

x = np.arange(100, 250)/100
energy = []
separation = []
for a in x:
    mol.atom = '''La 0 0 0; H %.2f 0 0''' % a
    mol.build()
    mf = scf.RHF(mol).run()
    mycc = cc.CCSD(mf).run()
    r = np.copy(mol.atom_coords(unit='ANG'))
    separation.append(np.linalg.norm(r[0] - r[1]))
    energy.append(mycc.e_tot)
data = np.zeros((len(separation), 2))
data[:,0] = separation
data[:,1] = energy
np.savetxt('CCSD.dat', data, fmt='%.5f')

where path/to/def2-qzvppd.1.nw should be the La ECP. This ECP comes from Basis Set Exchange, and I will update the repo with a link to it. The energy is written to CCSD.dat, which I inspected and found satisfactory. I then passed it to this script:

import numpy as np
import matplotlib.pyplot as plt
import sys
from ase.units import Hartree

# name of energy data
# has two columns: dimer separation, energy
fname = sys.argv[1]

data = np.loadtxt(fname)
data[:,1] *= Hartree # convert to eV
data[:,1] -= np.min(data[:,1]) # shift so that the minimum is at zero
cutoff = np.argmin(data[:,1]) # only keep the data up to the minimum
force = -np.gradient(data[:,1], data[:,0])

# plot just for sanity check
print('plotting data...')
plt.plot(data[1:cutoff,0], data[1:cutoff,1], 'b.-')
plt.xlabel('r (angstrom)')
plt.ylabel('energy (eV)')
plt.grid()
plt.twinx()
plt.plot(data[1:cutoff,0], force[1:cutoff], 'g.-')
plt.ylabel('force (eV / angstrom)')
plt.title('energy (blue) and force (green)')
plt.tight_layout()
plt.show()

# write table
N = cutoff - 1
lines = []
lines.append('# generated from ' + fname + '\n')
lines.append('\n')
lines.append('your keyword here\n')
lines.append('N %s\n' % N)
lines.append('\n')
for i in range(1, cutoff):
    lines.append('{0:<3} {1:.6f} {2:.6f} {3:.6f}\n'.format(i, data[i,0],\
            data[i,1], force[i]))
with open('test.table', 'w') as f:
    f.writelines(lines)

which writes the final table to test.table. You'll notice that, in addition to converting the energy to eV, I have also shifted it up so that it goes to 0 at the end.

The La-La and H-H energies can be computed a similar way.

kkly1995 / Data-for-fcc-LaH

pyscf #1