Skip to content

re-enable gaussian MM charges #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fishjojo merged 7 commits into from
Sep 27, 2025
Merged

re-enable gaussian MM charges #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
27a8af5
re-enable gaussian MM charges
MoleOrbitalHybridAnalyst Sep 25, 2025
8ae05a8
add periodic QMMM using multipoles
MoleOrbitalHybridAnalyst Sep 25, 2025
7dd4f16
minor
MoleOrbitalHybridAnalyst Sep 25, 2025
711704a
minor
MoleOrbitalHybridAnalyst Sep 25, 2025
bd6d4df
pep8
MoleOrbitalHybridAnalyst Sep 25, 2025
39f64f1
pep8
MoleOrbitalHybridAnalyst Sep 25, 2025
5d085b7
pep8
MoleOrbitalHybridAnalyst Sep 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 26 additions & 5 deletions pyscf/qmmm/itrf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@
from pyscf.lib import logger
from pyscf.qmmm import mm_mole

from scipy.special import erf


def add_mm_charges(scf_method, atoms_or_coords, charges, radii=None, unit=None):
'''Embedding the one-electron (non-relativistic) potential generated by MM
Expand Down Expand Up @@ -179,14 +181,18 @@ def get_hcore(self, mol=None):

def energy_nuc(self):
# interactions between QM nuclei and MM particles
assert self.mm_mol.charge_model != 'gaussian' # TODO: support Gaussian charge, same as the one in PCM
nuc = self.mol.energy_nuc()
coords = self.mm_mol.atom_coords()
charges = self.mm_mol.atom_charges()
if self.mm_mol.charge_model == 'gaussian':
expnts = numpy.sqrt(self.mm_mol.get_zetas())
for j in range(self.mol.natm):
q2, r2 = self.mol.atom_charge(j), self.mol.atom_coord(j)
r = lib.norm(r2-coords, axis=1)
nuc += q2*(charges/r).sum()
if self.mm_mol.charge_model != 'gaussian':
nuc += q2*(charges/r).sum()
else:
nuc += q2*(charges*erf(expnts*r)/r).sum()
return nuc

def to_gpu(self):
Expand Down Expand Up @@ -371,9 +377,11 @@ def grad_hcore_mm(self, dm, mol=None):

def grad_nuc(self, mol=None, atmlst=None):
if mol is None: mol = self.mol
assert self.base.mm_mol.charge_model != 'gaussian' # TODO: support Gaussian charge, same as the one in PCM
coords = self.base.mm_mol.atom_coords()
charges = self.base.mm_mol.atom_charges()
if self.base.mm_mol.charge_model == 'gaussian':
expnts = self.base.mm_mol.get_zetas()
radii = 1 / numpy.sqrt(expnts)

g_qm = super().grad_nuc(mol, atmlst)
# nuclei lattice interaction
Expand All @@ -382,7 +390,12 @@ def grad_nuc(self, mol=None, atmlst=None):
q1 = mol.atom_charge(i)
r1 = mol.atom_coord(i)
r = lib.norm(r1-coords, axis=1)
g_mm[i] = -q1 * numpy.einsum('i,ix,i->x', charges, r1-coords, 1/r**3)
if self.base.mm_mol.charge_model != 'gaussian':
coulkern = 1/r**3
else:
coulkern = erf(r/radii)/r - 2/(numpy.sqrt(numpy.pi)*radii) * numpy.exp(-expnts*r**2)
coulkern = coulkern / r**2
g_mm[i] = -q1 * numpy.einsum('i,ix,i->x', charges, r1-coords, coulkern)
if atmlst is not None:
g_mm = g_mm[atmlst]
return g_qm + g_mm
Expand All @@ -397,12 +410,20 @@ def grad_nuc_mm(self, mol=None):
mm_mol = self.base.mm_mol
coords = mm_mol.atom_coords()
charges = mm_mol.atom_charges()
if mm_mol.charge_model == 'gaussian':
expnts = mm_mol.get_zetas()
radii = 1 / numpy.sqrt(expnts)
g_mm = numpy.zeros_like(coords)
for i in range(mol.natm):
q1 = mol.atom_charge(i)
r1 = mol.atom_coord(i)
r = lib.norm(r1-coords, axis=1)
g_mm += q1 * numpy.einsum('i,ix,i->ix', charges, r1-coords, 1/r**3)
if mm_mol.charge_model != 'gaussian':
coulkern = 1/r**3
else:
coulkern = erf(r/radii)/r - 2/(numpy.sqrt(numpy.pi)*radii) * numpy.exp(-expnts*r**2)
coulkern = coulkern / r**2
g_mm += q1 * numpy.einsum('i,ix,i->ix', charges, r1-coords, coulkern)
return g_mm

def to_gpu(self):
Expand Down
1 change: 1 addition & 0 deletions pyscf/qmmm/pbc/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
from pyscf.qmmm.pbc import mm_mole, itrf
Loading