pyscf

[13yong]

Dear author:
Hello the author, I am a new man to study pyscf. Your article mentions using pyscf to generate the pair potential, but I followed the parameters provided in your article and was not able to reproduce your pair potential data. If it is convenient for you, can you post the script of pyscf, or upload it to github?
Looking forward to your reply

[kkly1995]

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.