34
34
35
35
from cupyx .scipy .special import erfc , erf
36
36
37
- def add_mm_charges (scf_method , atoms_or_coords , a , charges , radii = None ,
37
+ def add_mm_charges (scf_method , atoms_or_coords , a , charges , radii = None ,
38
38
rcut_ewald = None , rcut_hcore = None , unit = None ):
39
39
'''Embedding the one-electron (non-relativistic) potential generated by MM
40
40
point charges into QM Hamiltonian.
@@ -43,8 +43,8 @@ def add_mm_charges(scf_method, atoms_or_coords, a, charges, radii=None,
43
43
the nuclei in QM region and the MM charges, and the static Coulomb
44
44
interaction between the electron density and the MM charges. The electrostatic
45
45
interactions between reference cell and periodic images are also computed at
46
- point charge level. It does not include the static Coulomb interactions
47
- of the MM point charges, the MM energy, the vdw interaction or other
46
+ point charge level. It does not include the static Coulomb interactions
47
+ of the MM point charges, the MM energy, the vdw interaction or other
48
48
bonding/non-bonding effects between QM region and MM particles.
49
49
50
50
Args:
@@ -79,7 +79,7 @@ def add_mm_charges(scf_method, atoms_or_coords, a, charges, radii=None,
79
79
mol = scf_method .mol
80
80
if unit is None :
81
81
unit = mol .unit
82
- mm_mol = mm_mole .create_mm_mol (atoms_or_coords , a , charges , radii = radii ,
82
+ mm_mol = mm_mole .create_mm_mol (atoms_or_coords , a , charges , radii = radii ,
83
83
rcut_ewald = rcut_ewald , rcut_hcore = rcut_hcore , unit = unit )
84
84
return qmmm_for_scf (scf_method , mm_mol )
85
85
@@ -117,7 +117,8 @@ class QMMM:
117
117
class QMMMSCF (QMMM ):
118
118
_keys = {'mm_mol' , 's1r' , 's1rr' , 'mm_ewald_pot' , 'qm_ewald_hess' , 'e_nuc' }
119
119
120
- to_cpu = NotImplemented
120
+ to_cpu = NotImplemented
121
+ as_scanner = NotImplemented
121
122
122
123
def __init__ (self , method , mm_mol ):
123
124
self .__dict__ .update (method .__dict__ )
@@ -304,7 +305,7 @@ def get_qm_quadrupoles(self, dm, s1rr=None):
304
305
def get_vdiff (self , mol , ewald_pot ):
305
306
'''
306
307
vdiff_uv = d Q_I / d dm_uv ewald_pot[0]_I
307
- + d D_Ix / d dm_uv ewald_pot[1]_Ix
308
+ + d D_Ix / d dm_uv ewald_pot[1]_Ix
308
309
+ d O_Ixy / d dm_uv ewald_pot[2]_Ixy
309
310
'''
310
311
vdiff = cp .zeros ((mol .nao , mol .nao ))
@@ -406,7 +407,7 @@ def energy_nuc(self):
406
407
mol = self .mol
407
408
Ls = self .mm_mol .get_lattice_Ls ()
408
409
qm_center = np .mean (mol .atom_coords (), axis = 0 )
409
- all_coords = lib .direct_sum ('ix+Lx->Lix' ,
410
+ all_coords = lib .direct_sum ('ix+Lx->Lix' ,
410
411
self .mm_mol .atom_coords (), Ls ).reshape (- 1 ,3 )
411
412
all_charges = np .hstack ([self .mm_mol .atom_charges ()] * len (Ls ))
412
413
all_expnts = np .hstack ([np .sqrt (self .mm_mol .get_zetas ())] * len (Ls ))
@@ -441,7 +442,7 @@ def nuc_grad_method(self):
441
442
Gradients = nuc_grad_method
442
443
443
444
444
- def add_mm_charges_grad (scf_grad , atoms_or_coords , a , charges , radii = None ,
445
+ def add_mm_charges_grad (scf_grad , atoms_or_coords , a , charges , radii = None ,
445
446
rcut_ewald = None , rcut_hcore = None , unit = None ):
446
447
'''Apply the MM charges in the QM gradients' method. It affects both the
447
448
electronic and nuclear parts of the QM fragment.
@@ -473,7 +474,7 @@ def add_mm_charges_grad(scf_grad, atoms_or_coords, a, charges, radii=None,
473
474
mol = scf_grad .mol
474
475
if unit is None :
475
476
unit = mol .unit
476
- mm_mol = mm_mole .create_mm_mol (atoms_or_coords , a , charges , radii = radii ,
477
+ mm_mol = mm_mole .create_mm_mol (atoms_or_coords , a , charges , radii = radii ,
477
478
rcut_ewald = rcut_ewald , rcut_hcore = rcut_hcore , unit = unit )
478
479
mm_grad = qmmm_grad_for_scf (scf_grad )
479
480
mm_grad .base .mm_mol = mm_mol
@@ -581,7 +582,7 @@ def grad_ewald(self, dm=None, with_mm=False, mm_ewald_pot=None, qm_ewald_pot=Non
581
582
s1r .append (
582
583
cp .asarray (- mol .intor ('int1e_irp' , shls_slice = shlslc ).
583
584
reshape (3 , 3 , - 1 , mol .nao )))
584
- # s1rr[a,b,x,u,v] =
585
+ # s1rr[a,b,x,u,v] =
585
586
# \int phi_u [3/2*(r_a-Ri_a)(r_b-Ri_b)-1/2*(r-Ri)^2 delta_ab] (-\nable_x phi_v) dr
586
587
s1rr_ = cp .asarray (- mol .intor ('int1e_irrp' , shls_slice = shlslc ).
587
588
reshape (3 , 3 , 3 , - 1 , mol .nao ))
@@ -635,8 +636,8 @@ def grad_ewald(self, dm=None, with_mm=False, mm_ewald_pot=None, qm_ewald_pot=Non
635
636
grad_Tij = lambda R , r : get_multipole_tensors_pp (R , [1 ,2 ,3 ], r )
636
637
grad_kTij = lambda R , r , eta : get_multipole_tensors_pg (R , eta , [1 ,2 ,3 ], r )
637
638
638
- def grad_qm_multipole (Tija , Tijab , Tijabc ,
639
- qm_charges , qm_dipoles , qm_quads ,
639
+ def grad_qm_multipole (Tija , Tijab , Tijabc ,
640
+ qm_charges , qm_dipoles , qm_quads ,
640
641
mm_charges ):
641
642
Tc = contract ('ijx,j->ix' , Tija , mm_charges )
642
643
res = contract ('i,ix->ix' , qm_charges , Tc )
@@ -646,8 +647,8 @@ def grad_qm_multipole(Tija, Tijab, Tijabc,
646
647
res += contract ('iab,ixab->ix' , qm_quads , Tc ) / 3
647
648
return res
648
649
649
- def grad_mm_multipole (Tija , Tijab , Tijabc ,
650
- qm_charges , qm_dipoles , qm_quads ,
650
+ def grad_mm_multipole (Tija , Tijab , Tijabc ,
651
+ qm_charges , qm_dipoles , qm_quads ,
651
652
mm_charges ):
652
653
Tc = contract ('i,ijx->jx' , qm_charges , Tija )
653
654
Tc += contract ('ia,ijxa->jx' , qm_dipoles , Tijab )
@@ -771,7 +772,7 @@ def grad_mm_multipole(Tija, Tijab, Tijabc,
771
772
# ---------------------------------------------- #
772
773
# ---------- Ewald k-space gradient ------------ #
773
774
# ---------------------------------------------- #
774
-
775
+
775
776
mesh = cell .mesh
776
777
Gv , Gvbase , weights = cell .get_Gv_weights (mesh )
777
778
Gv = cp .asarray (Gv )
@@ -967,7 +968,7 @@ def grad_mm_multipole(Tija, Tijab, Tijabc,
967
968
temp = contract ('ag,bgx->abgx' , temp , temp2 )
968
969
temp = contract ('jg,abgx->abjx' , cosGvRqm , temp )
969
970
qm_ewg_grad -= contract ('jab,abjx->jx' , qm_quads , temp ) / 3
970
-
971
+
971
972
logger .timer (self , 'grad_ewald k-space' , * cput2 )
972
973
logger .timer (self , 'grad_ewald' , * cput0 )
973
974
if not with_mm :
@@ -1113,7 +1114,7 @@ def grad_nuc_mm(self, mol=None):
1113
1114
charges = mm_mol .atom_charges ()
1114
1115
Ls = mm_mol .get_lattice_Ls ()
1115
1116
qm_center = np .mean (mol .atom_coords (), axis = 0 )
1116
- all_coords = lib .direct_sum ('ix+Lx->Lix' ,
1117
+ all_coords = lib .direct_sum ('ix+Lx->Lix' ,
1117
1118
mm_mol .atom_coords (), Ls ).reshape (- 1 ,3 )
1118
1119
all_charges = np .hstack ([mm_mol .atom_charges ()] * len (Ls ))
1119
1120
all_expnts = np .hstack ([np .sqrt (mm_mol .get_zetas ())] * len (Ls ))
@@ -1130,7 +1131,7 @@ def grad_nuc_mm(self, mol=None):
1130
1131
r1 = mol .atom_coord (i )
1131
1132
r = lib .norm (r1 - coords , axis = 1 )
1132
1133
g_mm [mask ] += q1 * lib .einsum ('i,ix,i->ix' , charges , r1 - coords , erf (expnts * r )/ r ** 3 )
1133
- g_mm [mask ] -= q1 * lib .einsum ('i,ix,i->ix' , charges * expnts * 2 / np .sqrt (np .pi ),
1134
+ g_mm [mask ] -= q1 * lib .einsum ('i,ix,i->ix' , charges * expnts * 2 / np .sqrt (np .pi ),
1134
1135
r1 - coords , np .exp (- expnts ** 2 * r ** 2 )/ r ** 2 )
1135
1136
g_mm = g_mm .reshape (len (Ls ), - 1 , 3 )
1136
1137
g_mm = np .sum (g_mm , axis = 0 )
0 commit comments