diff --git a/dmff/admp/qeq.py b/dmff/admp/qeq.py
new file mode 100644
index 000000000..0c3e0cf77
--- /dev/null
+++ b/dmff/admp/qeq.py
@@ -0,0 +1,213 @@
+#!/usr/bin/env python
+import sys
+import absl 
+import numpy as np 
+import jax.numpy as jnp
+import openmm.app as app
+import openmm.unit as unit
+from dmff.settings import DO_JIT
+from dmff.common.constants import DIELECTRIC
+from dmff.common import nblist 
+from jax_md import space, partition
+from jax import grad, value_and_grad, vmap, jit
+from jaxopt import OptaxSolver
+from itertools import combinations
+import jaxopt
+import jax
+import scipy
+import pickle
+
+from jax.scipy.special import erf, erfc
+
+from dmff.utils import jit_condition, regularize_pairs, pair_buffer_scales
+
+
+jax.config.update("jax_enable_x64", True)
+
+class ADMPQeqForce:
+
+    def __init__(self, q, lagmt, damp_mod=3, neutral_flag=True, slab_flag=False,  constQ=True, pbc_flag = True):
+        
+        self.damp_mod = damp_mod
+        self.neutral_flag = neutral_flag
+        self.slab_flag = slab_flag
+        self.constQ = constQ
+        self.pbc_flag = pbc_flag
+        self.q = q
+        self.lagmt = lagmt
+        return
+
+    def generate_get_energy(self):
+      #  q = self.q
+        damp_mod = self.damp_mod
+        neutral_flag = self.neutral_flag
+        constQ = self.constQ
+        pbc_flag = self.pbc_flag
+      #  lagmt = self.lagmt
+
+        if eval(constQ) is True:
+            e_constraint = E_constQ
+        else:
+            e_constraint = E_constP
+        self.e_constraint = e_constraint
+
+        if eval(damp_mod) is False:
+            e_sr = E_sr0
+            e_site = E_site
+        elif eval(damp_mod) == 2:
+            e_sr = E_sr2
+            e_site = E_site2
+        elif eval(damp_mod) == 3:
+            e_sr = E_sr3
+            e_site = E_site3
+
+      #  if pbc_flag is False:
+      #      e_coul = E_CoulNocutoff
+      #  else:
+      #      e_coul = E_coul 
+        def get_energy(positions, box, pairs, q, lagmt, eta, chi, J, const_list, const_vals,pme_generator):
+            
+            pos = positions
+            ds = ds_pairs(pos, box, pairs, pbc_flag)
+            buffer_scales = pair_buffer_scales(pairs)
+            kappa = pme_generator.coulforce.kappa
+            def E_full(q, lagmt, const_vals, chi, J, pos, box, pairs, eta, ds, buffer_scales):
+                e1 = e_constraint(q, lagmt, const_list, const_vals)
+                e2 = e_sr(pos*10, box*10 ,pairs , q , eta, ds*10, buffer_scales)
+                e3 = e_site( chi, J , q)
+                e4 = pme_generator.coulenergy(pos, box ,pairs, q, pme_generator.mscales_coul)
+                e5 = E_corr(pos*10, box*10, pairs, q, kappa/10, neutral_flag)
+                return e1 + e2 + e3 + e4 + e5 
+            @jit
+            def E_grads(b_value, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales):
+                n_const = len(const_vals)
+                q = b_value[:-n_const]
+                lagmt = b_value[-n_const:]
+                g1,g2 = grad(E_full,argnums=(0,1))(q, lagmt, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales)
+                g = jnp.concatenate((g1,g2))
+                return g
+            
+            def Q_equi(b_value, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales):
+                rf=jaxopt.ScipyRootFinding(optimality_fun=E_grads,method='hybr',jit=False,tol=1e-10)
+                q0,state1 = rf.run(b_value, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales)
+                return q0,state1
+    
+            def get_chgs():
+                n_const = len(self.lagmt)
+                b_value = jnp.concatenate((self.q,self.lagmt))
+                q0,state1 = Q_equi(b_value, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales)
+                self.q = q0[:-n_const]
+                self.lagmt = q0[-n_const:]
+                return q0,state1
+
+            q0,state1 = get_chgs()
+            self.q0 = q0
+            self.state1 = state1
+            energy = E_full(self.q, self.lagmt, const_vals, chi, J, positions, box, pairs, eta, ds , buffer_scales)
+            self.e_grads = E_grads(q0, const_vals, chi, J, positions, box, pairs, eta, ds, buffer_scales)
+            self.e_full = E_full
+            return  energy 
+
+        return get_energy
+    def update_env(self, attr, val):
+        '''
+        Update the environment of the calculator
+        '''
+        setattr(self, attr, val)
+        self.refresh_calculators()
+
+
+    def refresh_calculators(self):
+        '''
+        refresh the energy and force calculators according to the current environment
+        '''
+        # generate the force calculator
+        self.get_energy = self.generate_get_energy()
+        self.get_forces = value_and_grad(self.get_energy)
+        return
+    
+def E_constQ(q, lagmt, const_list, const_vals):
+    constraint = (jnp.sum(q[const_list], axis=1) - const_vals) * lagmt
+    return np.sum(constraint)
+def E_constP(q, lagmt, const_list, const_vals):
+    constraint = jnp.sum(q[const_list], axis=1) * const_vals
+    return np.sum(constraint)
+
+def E_sr(pos, box, pairs, q, eta, ds, buffer_scales ):
+    return 0
+def E_sr2(pos, box, pairs, q, eta, ds, buffer_scales ):
+    etasqrt = jnp.sqrt( 2 * ( jnp.array(eta)[pairs[:,0]] **2 + jnp.array(eta)[pairs[:,1]] **2))
+    pre_pair = - eta_piecewise(etasqrt,ds) * DIELECTRIC
+    pre_self = etainv_piecewise(eta) /( jnp.sqrt(2 * jnp.pi)) * DIELECTRIC
+    e_sr_pair = pre_pair * q[pairs[:,0]] * q[pairs[:,1]] /ds * buffer_scales
+    e_sr_self = pre_self * q * q
+    e_sr = jnp.sum(e_sr_pair) + jnp.sum(e_sr_self)
+    return e_sr
+def E_sr3(pos, box, pairs, q, eta, ds, buffer_scales ):
+    etasqrt = jnp.sqrt( jnp.array(eta)[pairs[:,0]] **2 +  jnp.array(eta)[pairs[:,1]] **2 )
+    pre_pair = - eta_piecewise(etasqrt,ds) * DIELECTRIC
+    pre_self = etainv_piecewise(eta) /( jnp.sqrt(2 * jnp.pi)) * DIELECTRIC
+    e_sr_pair = pre_pair * q[pairs[:,0]] * q[pairs[:,1]] /ds * buffer_scales
+    e_sr_self = pre_self * q * q
+    e_sr = jnp.sum(e_sr_pair) + jnp.sum(e_sr_self)
+    return e_sr
+
+def E_site(chi, J , q ):
+    return 0
+def E_site2(chi, J , q ):
+    ene = (chi * q + 0.5 * J * q **2 ) * 96.4869
+    return np.sum(ene)
+def E_site3(chi, J , q ):
+    ene =  chi * q *4.184 + J * q **2 *DIELECTRIC * 2 * jnp.pi
+    return np.sum(ene)
+
+def E_corr(pos, box, pairs, q, kappa, neutral_flag = True):
+   # def E_corr():
+    V = jnp.linalg.det(box)
+    pre_corr = 2 * jnp.pi / V * DIELECTRIC
+    Mz = jnp.sum(q * pos[:,2])
+    Q_tot = jnp.sum(q)
+    Lz = jnp.linalg.norm(box[3])
+    e_corr = pre_corr * (Mz **2 - Q_tot * (jnp.sum(q * pos[:,2] **2)) - Q_tot **2 * Lz **2 /12)
+    if eval(neutral_flag) is True:
+      #  kappa = pme_potential.pme_force.kappa
+        pre_corr_non = - jnp.pi / (2 * V * kappa **2) * DIELECTRIC
+        e_corr_non = pre_corr_non * Q_tot **2 
+        e_corr += e_corr_non
+    return np.sum( e_corr)
+
+def E_CoulNocutoff(pos, box, pairs, q, ds):
+    e = q[pairs[:,0]] * q[pairs[:,1]] /ds * DIELECTRIC
+    return jnp.sum(e)
+
+def E_Coul(pos, box, pairs, q, ds):
+    return 0 
+
+@jit_condition(static_argnums=(3))
+def ds_pairs(positions, box, pairs, pbc_flag):
+    pos1 = positions[pairs[:,0].astype(int)]
+    pos2 = positions[pairs[:,1].astype(int)]
+    if pbc_flag is False:
+        dr = pos1 - pos2
+    else:
+        box_inv = jnp.linalg.inv(box)
+        dpos = pos1 - pos2
+        dpos = dpos.dot(box_inv)
+        dpos -= jnp.floor(dpos+0.5)
+        dr = dpos.dot(box)
+    ds = jnp.linalg.norm(dr,axis=1)
+    return ds
+
+@jit_condition()
+@vmap
+def eta_piecewise(eta,ds):
+    return jnp.piecewise(eta, (eta > 1e-4, eta <= 1e-4),
+                        (lambda x: jnp.array(erfc( ds / eta)), lambda x:jnp.array(0))) 
+                         
+@jit_condition()
+@vmap
+def etainv_piecewise(eta):
+    return jnp.piecewise(eta, (eta > 1e-4, eta <= 1e-4),
+                        (lambda x: jnp.array(1/eta), lambda x:jnp.array(0))) 
+                         
+
diff --git a/dmff/generators/QeqGenerator.py b/dmff/generators/QeqGenerator.py
new file mode 100644
index 000000000..aba48d0f8
--- /dev/null
+++ b/dmff/generators/QeqGenerator.py
@@ -0,0 +1,135 @@
+#!/usr/bin/env python
+
+import openmm.app as app
+import openmm.unit as unit
+from typing import Tuple
+import numpy as np
+import jax.numpy as jnp
+import jax
+from dmff.api.topology import DMFFTopology
+from dmff.api.paramset import ParamSet
+from dmff.api.xmlio import XMLIO
+from dmff.api.hamiltonian import _DMFFGenerators
+from dmff.utils import DMFFException, isinstance_jnp
+from dmff.admp.qeq import ADMPQeqForce 
+from dmff.generators.classical import CoulombGenerator
+from dmff.admp import qeq 
+
+
+class ADMPQeqGenerator:
+    def __init__(self, ffinfo:dict, paramset: ParamSet):
+
+        self.name = 'ADMPQeqForce'
+        self.ffinfo = ffinfo
+        paramset.addField(self.name)
+        self.key_type = None
+        keys , params = [], []
+        for node in self.ffinfo["Forces"][self.name]["node"]:
+            attribs = node["attrib"]
+            
+            if self.key_type is None and "type" in attribs:
+                self.key_type = "type"
+            elif self.key_type is None and "class" in attribs:
+                self.key_type = "class"
+            elif self.key_type is not None and f"{self.key_type}" not in attribs:
+                raise ValueError("Keyword 'class' or 'type' cannot be used together.")
+            elif self.key_type is not None and f"{self.key_type}" in attribs:
+                pass
+            else:
+                raise ValueError("Cannot find key type for ADMPQeqForce.")
+            key = attribs[self.key_type]
+            keys.append(key)
+
+            chi0 = float(attribs["chi"]) 
+            J0 = float(attribs["J"])
+            eta0 = float(attribs["eta"])
+
+            params.append([chi0, J0, eta0])
+
+        self.keys = keys
+        chi = jnp.array([i[0] for i in params])
+        J = jnp.array([i[1] for i in params])
+        eta = jnp.array([i[2] for i in params])
+
+        paramset.addParameter(chi, "chi", field=self.name)
+        paramset.addParameter(J, "J", field=self.name)
+        paramset.addParameter(eta, "eta", field=self.name)
+        # default params
+        self._jaxPotential = None
+        self.damp_mod = self.ffinfo["Forces"][self.name]["meta"]["DampMod"] 
+        self.neutral_flag = self.ffinfo["Forces"][self.name]["meta"]["NeutralFlag"]
+        self.slab_flag = self.ffinfo["Forces"][self.name]["meta"]["SlabFlag"]
+        self.constQ = self.ffinfo["Forces"][self.name]["meta"]["ConstQFlag"]
+        self.pbc_flag = self.ffinfo["Forces"][self.name]["meta"]["PbcFlag"] 
+
+        self.pme_generator = CoulombGenerator(ffinfo, paramset)
+
+    def getName(self) -> str:
+        """
+        Returns the name of the force field.
+
+        Returns:
+        --------
+        str
+            The name of the force field.
+        """
+        return self.name
+    
+    def overwrite(self, paramset:ParamSet) -> None:
+
+        node_indices = [ i for i in range(len(self.ffinfo["Forces"][self.name]["node"])) if self.ffinfo["Forces"][self.name]["node"][i]["name"] == "QeqAtom"]
+        chi = paramset[self.name]["chi"]
+        J = paramset[self.name]["J"]
+        eta = paramset[self.name]["eta"]
+        for nnode, key in enumerate(self.keys):
+            self.ffinfo["Forces"][self.name]["node"][node_indices[nnode]]["attrib"] = {}
+            self.ffinfo["Forces"][self.name]["node"][node_indices[nnode]]["attrib"][f"{self.key_type}"] = key
+            chi0 = chi[nnode]
+            J0 = J[nnode]
+            eta0 = eta[nnode]
+            self.ffinfo["Forces"][self.name]["node"][bond_node_indices[nnode]]["attrib"]["chi"] = str(chi0)
+            self.ffinfo["Forces"][self.name]["node"][bond_node_indices[nnode]]["attrib"]["J"] = str(J0)
+            self.ffinfo["Forces"][self.name]["node"][bond_node_indices[nnode]]["attrib"]["eta"] = str(eta0)
+
+
+    def createPotential(self, topdata: DMFFTopology, nonbondedMethod, nonbondedCutoff, charges, const_list, const_vals, map_atomtype):
+
+        n_atoms = topdata._numAtoms
+        n_residues = topdata._numResidues
+        
+        q = jnp.array(charges)
+        lagmt = np.ones(n_residues)
+        b_value = jnp.concatenate((q,lagmt))
+        qeq_force = ADMPQeqForce(q, lagmt,self.damp_mod, self.neutral_flag,
+                                 self.slab_flag, self.constQ, self.pbc_flag)
+        self.qeq_force = qeq_force
+        qeq_energy = qeq_force.generate_get_energy()
+
+        self.pme_potential = self.pme_generator.createPotential(topdata, app.PME, nonbondedCutoff )  
+        def potential_fn(positions: jnp.ndarray, box: jnp.ndarray, pairs: jnp.ndarray, params: ParamSet) -> jnp.ndarray:
+            
+            n_atoms = len(positions)
+           # map_atomtype = np.zeros(n_atoms)
+            eta = np.array(params[self.name]["eta"])[map_atomtype]
+            chi = np.array(params[self.name]["chi"])[map_atomtype]
+            J = np.array(params[self.name]["J"])[map_atomtype]
+            self.eta = jnp.array(eta)
+            self.chi = jnp.array(chi)
+            self.J = jnp.array(J)
+         #   coulenergy = self.pme_generator.coulenergy
+         #   pme_energy = pme_potential(positions, box, pairs, params)
+            damp_mod = self.damp_mod
+            neutral_flag = self.neutral_flag
+            constQ = self.constQ
+            pbc_flag = self.pbc_flag
+            
+            qeq_energy0 = qeq_energy(positions, box, pairs, q, lagmt,
+                                    eta, chi, J,const_list, 
+                                    const_vals, self.pme_generator) 
+           # return pme_energy + qeq_energy0
+            return qeq_energy0
+
+        self._jaxPotential = potential_fn
+        return potential_fn
+
+_DMFFGenerators["ADMPQeqForce"] = ADMPQeqGenerator
diff --git a/dmff/generators/classical.py b/dmff/generators/classical.py
index 620d0a738..e6456d94f 100644
--- a/dmff/generators/classical.py
+++ b/dmff/generators/classical.py
@@ -1054,6 +1054,7 @@ def createPotential(self, topdata: DMFFTopology, nonbondedMethod,
         mscales_coul = jnp.array([0.0, 0.0, 0.0, 1.0, 1.0,
                                   1.0])  # mscale for PME
         mscales_coul = mscales_coul.at[2].set(self.coulomb14scale)
+        self.mscales_coul = mscales_coul # for qeq calculation
 
         # set PBC
         if nonbondedMethod not in [app.NoCutoff, app.CutoffNonPeriodic]:
@@ -1075,7 +1076,8 @@ def createPotential(self, topdata: DMFFTopology, nonbondedMethod,
         if nonbondedMethod is app.PME:
             cell = topdata.getPeriodicBoxVectors()
             box = jnp.array(cell)
-            self.ethresh = kwargs.get("ethresh", 1e-6)
+           # self.ethresh = kwargs.get("ethresh", 1e-6)
+            self.ethresh = kwargs.get("ethresh", 5e-4) #for qeq calculation
             self.coeff_method = kwargs.get("PmeCoeffMethod", "openmm")
             self.fourier_spacing = kwargs.get("PmeSpacing", 0.1)
             kappa, K1, K2, K3 = setup_ewald_parameters(r_cut, self.ethresh,
@@ -1120,7 +1122,8 @@ def createPotential(self, topdata: DMFFTopology, nonbondedMethod,
                 topology_matrix=top_mat if self._use_bcc else None)
 
         coulenergy = coulforce.generate_get_energy()
-
+        self.coulforce = coulforce  #for qeq calculation
+        self.coulenergy = coulenergy #for qeq calculation
         def potential_fn(positions, box, pairs, params):
 
             # check whether args passed into potential_fn are jnp.array and differentiable
diff --git a/tests/data/qeq.pdb b/tests/data/qeq.pdb
new file mode 100644
index 000000000..c835be211
--- /dev/null
+++ b/tests/data/qeq.pdb
@@ -0,0 +1,149 @@
+HEADER    electrode systems
+TITLE     MDANALYSIS FRAME 0: Created by PDBWriter
+REMARK     GENERATED BY Haichao
+CRYST1   22.116   17.184  200.000  90.00  90.00  90.00 P 1           1
+ATOM      1 C001 CNC X   11      2.460  14.913  10.100  1.00  0.00      C     
+ATOM      2 C002 CNC X   11      4.920  14.913  10.100  1.00  0.00      C     
+ATOM      3 C003 CNC X   11      7.380  14.913  10.100  1.00  0.00      C     
+ATOM      4 C004 CNC X   11      9.840  14.913  10.100  1.00  0.00      C     
+ATOM      5 C005 CNC X   11     12.300  14.913  10.100  1.00  0.00      C     
+ATOM      6 C006 CNC X   11     14.760  14.913  10.100  1.00  0.00      C     
+ATOM      7 C007 CNC X   11     17.220  14.913  10.100  1.00  0.00      C     
+ATOM      8 C008 CNC X   11      0.000  14.913  10.100  1.00  0.00      C     
+ATOM      9 C009 CNC X   11     19.680  14.913  10.100  1.00  0.00      C     
+ATOM     10 C010 CNC X   11     11.070   0.000  10.100  1.00  0.00      C     
+ATOM     11 C011 CNC X   11     12.300   0.710  10.100  1.00  0.00      C     
+ATOM     12 C012 CNC X   11     13.530   0.000  10.100  1.00  0.00      C     
+ATOM     13 C013 CNC X   11     14.760   0.710  10.100  1.00  0.00      C     
+ATOM     14 C014 CNC X   11     15.990   0.000  10.100  1.00  0.00      C     
+ATOM     15 C015 CNC X   11     17.220   0.710  10.100  1.00  0.00      C     
+ATOM     16 C016 CNC X   11     18.450   0.000  10.100  1.00  0.00      C     
+ATOM     17 C017 CNC X   11      0.000   0.710  10.100  1.00  0.00      C     
+ATOM     18 C018 CNC X   11     19.680   0.710  10.100  1.00  0.00      C     
+ATOM     19 C019 CNC X   11      1.230   0.000  10.100  1.00  0.00      C     
+ATOM     20 C020 CNC X   11     20.910   0.000  10.100  1.00  0.00      C     
+ATOM     21 C021 CNC X   11      2.460   0.710  10.100  1.00  0.00      C     
+ATOM     22 C022 CNC X   11      3.690   0.000  10.100  1.00  0.00      C     
+ATOM     23 C023 CNC X   11      4.920   0.710  10.100  1.00  0.00      C     
+ATOM     24 C024 CNC X   11      6.150   0.000  10.100  1.00  0.00      C     
+ATOM     25 C025 CNC X   11      7.380   0.710  10.100  1.00  0.00      C     
+ATOM     26 C026 CNC X   11      8.610   0.000  10.100  1.00  0.00      C     
+ATOM     27 C027 CNC X   11      9.840   0.710  10.100  1.00  0.00      C     
+ATOM     28 C028 CNC X   11      9.840   2.130  10.100  1.00  0.00      C     
+ATOM     29 C029 CNC X   11     11.070   2.841  10.100  1.00  0.00      C     
+ATOM     30 C030 CNC X   11     12.300   2.130  10.100  1.00  0.00      C     
+ATOM     31 C031 CNC X   11     13.530   2.841  10.100  1.00  0.00      C     
+ATOM     32 C032 CNC X   11     14.760   2.130  10.100  1.00  0.00      C     
+ATOM     33 C033 CNC X   11     15.990   2.841  10.100  1.00  0.00      C     
+ATOM     34 C034 CNC X   11     17.220   2.130  10.100  1.00  0.00      C     
+ATOM     35 C035 CNC X   11     18.450   2.841  10.100  1.00  0.00      C     
+ATOM     36 C036 CNC X   11      0.000   2.130  10.100  1.00  0.00      C     
+ATOM     37 C037 CNC X   11     19.680   2.130  10.100  1.00  0.00      C     
+ATOM     38 C038 CNC X   11      1.230   2.841  10.100  1.00  0.00      C     
+ATOM     39 C039 CNC X   11     20.910   2.841  10.100  1.00  0.00      C     
+ATOM     40 C040 CNC X   11      2.460   2.130  10.100  1.00  0.00      C     
+ATOM     41 C041 CNC X   11      3.690   2.841  10.100  1.00  0.00      C     
+ATOM     42 C042 CNC X   11      4.920   2.130  10.100  1.00  0.00      C     
+ATOM     43 C043 CNC X   11      6.150   2.841  10.100  1.00  0.00      C     
+ATOM     44 C044 CNC X   11      7.380   2.130  10.100  1.00  0.00      C     
+ATOM     45 C045 CNC X   11      8.610   2.841  10.100  1.00  0.00      C     
+ATOM     46 C046 CNC X   11      8.610   4.261  10.100  1.00  0.00      C     
+ATOM     47 C047 CNC X   11      9.840   4.971  10.100  1.00  0.00      C     
+ATOM     48 C048 CNC X   11     11.070   4.261  10.100  1.00  0.00      C     
+ATOM     49 C049 CNC X   11     12.300   4.971  10.100  1.00  0.00      C     
+ATOM     50 C050 CNC X   11     13.530   4.261  10.100  1.00  0.00      C     
+ATOM     51 C051 CNC X   11     14.760   4.971  10.100  1.00  0.00      C     
+ATOM     52 C052 CNC X   11     15.990   4.261  10.100  1.00  0.00      C     
+ATOM     53 C053 CNC X   11     17.220   4.971  10.100  1.00  0.00      C     
+ATOM     54 C054 CNC X   11     18.450   4.261  10.100  1.00  0.00      C     
+ATOM     55 C055 CNC X   11      0.000   4.971  10.100  1.00  0.00      C     
+ATOM     56 C056 CNC X   11     19.680   4.971  10.100  1.00  0.00      C     
+ATOM     57 C057 CNC X   11      1.230   4.261  10.100  1.00  0.00      C     
+ATOM     58 C058 CNC X   11     20.910   4.261  10.100  1.00  0.00      C     
+ATOM     59 C059 CNC X   11      2.460   4.971  10.100  1.00  0.00      C     
+ATOM     60 C060 CNC X   11      3.690   4.261  10.100  1.00  0.00      C     
+ATOM     61 C061 CNC X   11      4.920   4.971  10.100  1.00  0.00      C     
+ATOM     62 C062 CNC X   11      6.150   4.261  10.100  1.00  0.00      C     
+ATOM     63 C063 CNC X   11      7.380   4.971  10.100  1.00  0.00      C     
+ATOM     64 C064 CNC X   11      7.380   6.391  10.100  1.00  0.00      C     
+ATOM     65 C065 CNC X   11      8.610   7.101  10.100  1.00  0.00      C     
+ATOM     66 C066 CNC X   11      9.840   6.391  10.100  1.00  0.00      C     
+ATOM     67 C067 CNC X   11     12.300   6.391  10.100  1.00  0.00      C     
+ATOM     68 C068 CNC X   11     13.530   7.101  10.100  1.00  0.00      C     
+ATOM     69 C069 CNC X   11     14.760   6.391  10.100  1.00  0.00      C     
+ATOM     70 C070 CNC X   11     15.990   7.101  10.100  1.00  0.00      C     
+ATOM     71 C071 CNC X   11     17.220   6.391  10.100  1.00  0.00      C     
+ATOM     72 C072 CNC X   11     18.450   7.101  10.100  1.00  0.00      C     
+ATOM     73 C073 CNC X   11      0.000   6.391  10.100  1.00  0.00      C     
+ATOM     74 C074 CNC X   11     19.680   6.391  10.100  1.00  0.00      C     
+ATOM     75 C075 CNC X   11      1.230   7.101  10.100  1.00  0.00      C     
+ATOM     76 C076 CNC X   11     20.910   7.101  10.100  1.00  0.00      C     
+ATOM     77 C077 CNC X   11      2.460   6.391  10.100  1.00  0.00      C     
+ATOM     78 C078 CNC X   11      3.690   7.101  10.100  1.00  0.00      C     
+ATOM     79 C079 CNC X   11      4.920   6.391  10.100  1.00  0.00      C     
+ATOM     80 C080 CNC X   11      6.150   7.101  10.100  1.00  0.00      C     
+ATOM     81 C081 CNC X   11      6.150   8.522  10.100  1.00  0.00      C     
+ATOM     82 C082 CNC X   11      7.380   9.232  10.100  1.00  0.00      C     
+ATOM     83 C083 CNC X   11      8.610   8.522  10.100  1.00  0.00      C     
+ATOM     84 C084 CNC X   11      9.840   9.232  10.100  1.00  0.00      C     
+ATOM     85 C085 CNC X   11     11.070   8.522  10.100  1.00  0.00      C     
+ATOM     86 C086 CNC X   11     12.300   9.232  10.100  1.00  0.00      C     
+ATOM     87 C087 CNC X   11     13.530   8.522  10.100  1.00  0.00      C     
+ATOM     88 C088 CNC X   11     14.760   9.232  10.100  1.00  0.00      C     
+ATOM     89 C089 CNC X   11     15.990   8.522  10.100  1.00  0.00      C     
+ATOM     90 C090 CNC X   11     17.220   9.232  10.100  1.00  0.00      C     
+ATOM     91 C091 CNC X   11     18.450   8.522  10.100  1.00  0.00      C     
+ATOM     92 C092 CNC X   11      0.000   9.232  10.100  1.00  0.00      C     
+ATOM     93 C093 CNC X   11     19.680   9.232  10.100  1.00  0.00      C     
+ATOM     94 C094 CNC X   11      1.230   8.522  10.100  1.00  0.00      C     
+ATOM     95 C095 CNC X   11     20.910   8.522  10.100  1.00  0.00      C     
+ATOM     96 C096 CNC X   11      2.460   9.232  10.100  1.00  0.00      C     
+ATOM     97 C097 CNC X   11      3.690   8.522  10.100  1.00  0.00      C     
+ATOM     98 C098 CNC X   11      4.920   9.232  10.100  1.00  0.00      C     
+ATOM     99 C099 CNC X   11      4.920  10.652  10.100  1.00  0.00      C     
+ATOM    100 C100 CNC X   11      6.150  11.362  10.100  1.00  0.00      C     
+ATOM    101 C101 CNC X   11      7.380  10.652  10.100  1.00  0.00      C     
+ATOM    102 C102 CNC X   11      8.610  11.362  10.100  1.00  0.00      C     
+ATOM    103 C103 CNC X   11      9.840  10.652  10.100  1.00  0.00      C     
+ATOM    104 C104 CNC X   11     11.070  11.362  10.100  1.00  0.00      C     
+ATOM    105 C105 CNC X   11     12.300  10.652  10.100  1.00  0.00      C     
+ATOM    106 C106 CNC X   11     13.530  11.362  10.100  1.00  0.00      C     
+ATOM    107 C107 CNC X   11     14.760  10.652  10.100  1.00  0.00      C     
+ATOM    108 C108 CNC X   11     15.990  11.362  10.100  1.00  0.00      C     
+ATOM    109 C109 CNC X   11     17.220  10.652  10.100  1.00  0.00      C     
+ATOM    110 C110 CNC X   11     18.450  11.362  10.100  1.00  0.00      C     
+ATOM    111 C111 CNC X   11      0.000  10.652  10.100  1.00  0.00      C     
+ATOM    112 C112 CNC X   11     19.680  10.652  10.100  1.00  0.00      C     
+ATOM    113 C113 CNC X   11      1.230  11.362  10.100  1.00  0.00      C     
+ATOM    114 C114 CNC X   11     20.910  11.362  10.100  1.00  0.00      C     
+ATOM    115 C115 CNC X   11      2.460  10.652  10.100  1.00  0.00      C     
+ATOM    116 C116 CNC X   11      3.690  11.362  10.100  1.00  0.00      C     
+ATOM    117 C117 CNC X   11      3.690  12.783  10.100  1.00  0.00      C     
+ATOM    118 C118 CNC X   11      4.920  13.493  10.100  1.00  0.00      C     
+ATOM    119 C119 CNC X   11      6.150  12.783  10.100  1.00  0.00      C     
+ATOM    120 C120 CNC X   11      7.380  13.493  10.100  1.00  0.00      C     
+ATOM    121 C121 CNC X   11      8.610  12.783  10.100  1.00  0.00      C     
+ATOM    122 C122 CNC X   11      9.840  13.493  10.100  1.00  0.00      C     
+ATOM    123 C123 CNC X   11     11.070  12.783  10.100  1.00  0.00      C     
+ATOM    124 C124 CNC X   11     12.300  13.493  10.100  1.00  0.00      C     
+ATOM    125 C125 CNC X   11     13.530  12.783  10.100  1.00  0.00      C     
+ATOM    126 C126 CNC X   11     14.760  13.493  10.100  1.00  0.00      C     
+ATOM    127 C127 CNC X   11     15.990  12.783  10.100  1.00  0.00      C     
+ATOM    128 C128 CNC X   11     17.220  13.493  10.100  1.00  0.00      C     
+ATOM    129 C129 CNC X   11     18.450  12.783  10.100  1.00  0.00      C     
+ATOM    130 C130 CNC X   11      0.000  13.493  10.100  1.00  0.00      C     
+ATOM    131 C131 CNC X   11     19.680  13.493  10.100  1.00  0.00      C     
+ATOM    132 C132 CNC X   11      1.230  12.783  10.100  1.00  0.00      C     
+ATOM    133 C133 CNC X   11     20.910  12.783  10.100  1.00  0.00      C     
+ATOM    134 C134 CNC X   11      2.460  13.493  10.100  1.00  0.00      C     
+ATOM    135 C135 CNC X   11      1.230  15.623  10.100  1.00  0.00      C     
+ATOM    136 C136 CNC X   11      3.690  15.623  10.100  1.00  0.00      C     
+ATOM    137 C137 CNC X   11     20.910  15.623  10.100  1.00  0.00      C     
+ATOM    138 C138 CNC X   11      6.150  15.623  10.100  1.00  0.00      C     
+ATOM    139 C139 CNC X   11      8.610  15.623  10.100  1.00  0.00      C     
+ATOM    140 C140 CNC X   11     11.070  15.623  10.100  1.00  0.00      C     
+ATOM    141 C141 CNC X   11     13.530  15.623  10.100  1.00  0.00      C     
+ATOM    142 C142 CNC X   11     15.990  15.623  10.100  1.00  0.00      C     
+ATOM    143 C143 CNC X   11     18.450  15.623  10.100  1.00  0.00      C     
+ATOM    144 N001 CNC X   11     11.070   7.101  10.100  1.00  0.00      N     
+END
diff --git a/tests/data/qeq.xml b/tests/data/qeq.xml
new file mode 100644
index 000000000..664609905
--- /dev/null
+++ b/tests/data/qeq.xml
@@ -0,0 +1,197 @@
+<!-- Note that the atom types "name" are numeric - this is because of how AMOEBA/MPID assign anchor atoms -->
+<ForceField>
+ <AtomTypes>
+  <Type name="1" class="CW" element="C" mass="12.0107"/>
+  <Type name="2" class="NW" element="N" mass="14.0067"/>
+  <Type name="3" class="CW" element="C" mass="12.0107"/>
+  <Type name="4" class="CW" element="C" mass="12.0107"/>
+  <Type name="5" class="CW" element="C" mass="12.0107"/>
+ </AtomTypes>
+ <Residues>
+  <Residue name="CNC">
+   <Atom name="C001" type="1" charge="-0.001721"/>
+   <Atom name="C002" type="1" charge="0.001835"/>
+   <Atom name="C003" type="1" charge="-0.000278"/>
+   <Atom name="C004" type="1" charge="0.001665"/>
+   <Atom name="C005" type="1" charge="0.001717"/>
+   <Atom name="C006" type="1" charge="0.001210"/>
+   <Atom name="C007" type="1" charge="0.001117"/>
+   <Atom name="C008" type="1" charge="0.001007"/>
+   <Atom name="C009" type="1" charge="0.000622"/>
+   <Atom name="C010" type="1" charge="0.001496"/>
+   <Atom name="C011" type="1" charge="-0.001773"/>
+   <Atom name="C012" type="1" charge="0.001438"/>
+   <Atom name="C013" type="1" charge="-0.010743"/>
+   <Atom name="C014" type="1" charge="0.000499"/>
+   <Atom name="C015" type="1" charge="0.001085"/>
+   <Atom name="C016" type="1" charge="0.000817"/>
+   <Atom name="C017" type="1" charge="-0.007348"/>
+   <Atom name="C018" type="1" charge="-0.000574"/>
+   <Atom name="C019" type="1" charge="-0.001629"/>
+   <Atom name="C020" type="1" charge="0.000905"/>
+   <Atom name="C021" type="1" charge="0.000752"/>
+   <Atom name="C022" type="1" charge="-0.001020"/>
+   <Atom name="C023" type="1" charge="0.002377"/>
+   <Atom name="C024" type="1" charge="-0.002274"/>
+   <Atom name="C025" type="1" charge="-0.008115"/>
+   <Atom name="C026" type="1" charge="-0.000166"/>
+   <Atom name="C027" type="1" charge="-0.001606"/>
+   <Atom name="C028" type="1" charge="0.006095"/>
+   <Atom name="C029" type="1" charge="-0.014802"/>
+   <Atom name="C030" type="1" charge="0.003875"/>
+   <Atom name="C031" type="1" charge="0.000269"/>
+   <Atom name="C032" type="1" charge="0.002752"/>
+   <Atom name="C033" type="1" charge="-0.001540"/>
+   <Atom name="C034" type="1" charge="0.002408"/>
+   <Atom name="C035" type="1" charge="-0.009637"/>
+   <Atom name="C036" type="1" charge="0.004210"/>
+   <Atom name="C037" type="1" charge="0.001351"/>
+   <Atom name="C038" type="1" charge="-0.001683"/>
+   <Atom name="C039" type="1" charge="0.003371"/>
+   <Atom name="C040" type="1" charge="0.003562"/>
+   <Atom name="C041" type="1" charge="-0.011557"/>
+   <Atom name="C042" type="1" charge="0.004116"/>
+   <Atom name="C043" type="1" charge="-0.001749"/>
+   <Atom name="C044" type="1" charge="0.002514"/>
+   <Atom name="C045" type="1" charge="-0.001687"/>
+   <Atom name="C046" type="1" charge="0.011372"/>
+   <Atom name="C047" type="4" charge="-0.040468"/>
+   <Atom name="C048" type="5" charge="0.013502"/>
+   <Atom name="C049" type="4" charge="-0.040911"/>
+   <Atom name="C050" type="1" charge="0.011195"/>
+   <Atom name="C051" type="1" charge="-0.013165"/>
+   <Atom name="C052" type="1" charge="-0.000075"/>
+   <Atom name="C053" type="1" charge="0.002269"/>
+   <Atom name="C054" type="1" charge="0.001824"/>
+   <Atom name="C055" type="1" charge="-0.007915"/>
+   <Atom name="C056" type="1" charge="-0.001707"/>
+   <Atom name="C057" type="1" charge="0.000568"/>
+   <Atom name="C058" type="1" charge="0.004427"/>
+   <Atom name="C059" type="1" charge="0.001604"/>
+   <Atom name="C060" type="1" charge="0.000705"/>
+   <Atom name="C061" type="1" charge="0.003556"/>
+   <Atom name="C062" type="1" charge="-0.000351"/>
+   <Atom name="C063" type="1" charge="-0.013660"/>
+   <Atom name="C064" type="1" charge="0.006938"/>
+   <Atom name="C065" type="4" charge="-0.034317"/>
+   <Atom name="C066" type="3" charge="0.101458"/>
+   <Atom name="C067" type="3" charge="0.099507"/>
+   <Atom name="C068" type="4" charge="-0.034583"/>
+   <Atom name="C069" type="1" charge="0.007959"/>
+   <Atom name="C070" type="1" charge="-0.000108"/>
+   <Atom name="C071" type="1" charge="0.000794"/>
+   <Atom name="C072" type="1" charge="-0.011437"/>
+   <Atom name="C073" type="1" charge="0.003400"/>
+   <Atom name="C074" type="1" charge="0.001088"/>
+   <Atom name="C075" type="1" charge="-0.001834"/>
+   <Atom name="C076" type="1" charge="0.002119"/>
+   <Atom name="C077" type="1" charge="0.001460"/>
+   <Atom name="C078" type="1" charge="-0.011831"/>
+   <Atom name="C079" type="1" charge="0.001373"/>
+   <Atom name="C080" type="1" charge="0.001910"/>
+   <Atom name="C081" type="1" charge="0.001257"/>
+   <Atom name="C082" type="1" charge="-0.010480"/>
+   <Atom name="C083" type="5" charge="0.009487"/>
+   <Atom name="C084" type="4" charge="-0.032903"/>
+   <Atom name="C085" type="3" charge="0.096115"/>
+   <Atom name="C086" type="4" charge="-0.037532"/>
+   <Atom name="C087" type="5" charge="0.011060"/>
+   <Atom name="C088" type="1" charge="-0.012612"/>
+   <Atom name="C089" type="1" charge="0.002939"/>
+   <Atom name="C090" type="1" charge="0.000817"/>
+   <Atom name="C091" type="1" charge="0.002880"/>
+   <Atom name="C092" type="1" charge="-0.009572"/>
+   <Atom name="C093" type="1" charge="-0.001819"/>
+   <Atom name="C094" type="1" charge="0.002040"/>
+   <Atom name="C095" type="1" charge="0.004718"/>
+   <Atom name="C096" type="1" charge="0.002701"/>
+   <Atom name="C097" type="1" charge="0.002054"/>
+   <Atom name="C098" type="1" charge="0.000962"/>
+   <Atom name="C099" type="1" charge="0.001426"/>
+   <Atom name="C100" type="1" charge="-0.000006"/>
+   <Atom name="C101" type="1" charge="0.000386"/>
+   <Atom name="C102" type="1" charge="0.002711"/>
+   <Atom name="C103" type="1" charge="0.005549"/>
+   <Atom name="C104" type="1" charge="-0.015649"/>
+   <Atom name="C105" type="1" charge="0.007179"/>
+   <Atom name="C106" type="1" charge="0.003943"/>
+   <Atom name="C107" type="1" charge="0.001799"/>
+   <Atom name="C108" type="1" charge="-0.000509"/>
+   <Atom name="C109" type="1" charge="0.001458"/>
+   <Atom name="C110" type="1" charge="-0.010764"/>
+   <Atom name="C111" type="1" charge="0.003486"/>
+   <Atom name="C112" type="1" charge="0.001249"/>
+   <Atom name="C113" type="1" charge="0.001740"/>
+   <Atom name="C114" type="1" charge="0.003551"/>
+   <Atom name="C115" type="1" charge="-0.001178"/>
+   <Atom name="C116" type="1" charge="-0.012407"/>
+   <Atom name="C117" type="1" charge="0.003436"/>
+   <Atom name="C118" type="1" charge="0.000688"/>
+   <Atom name="C119" type="1" charge="0.002018"/>
+   <Atom name="C120" type="1" charge="-0.009659"/>
+   <Atom name="C121" type="1" charge="0.001063"/>
+   <Atom name="C122" type="1" charge="0.003147"/>
+   <Atom name="C123" type="1" charge="-0.000343"/>
+   <Atom name="C124" type="1" charge="0.002149"/>
+   <Atom name="C125" type="1" charge="-0.000510"/>
+   <Atom name="C126" type="1" charge="-0.008943"/>
+   <Atom name="C127" type="1" charge="-0.000310"/>
+   <Atom name="C128" type="1" charge="0.003428"/>
+   <Atom name="C129" type="1" charge="0.001463"/>
+   <Atom name="C130" type="1" charge="-0.004922"/>
+   <Atom name="C131" type="1" charge="0.000570"/>
+   <Atom name="C132" type="1" charge="-0.001153"/>
+   <Atom name="C133" type="1" charge="0.000384"/>
+   <Atom name="C134" type="1" charge="0.003119"/>
+   <Atom name="C135" type="1" charge="-0.000234"/>
+   <Atom name="C136" type="1" charge="-0.011229"/>
+   <Atom name="C137" type="1" charge="0.001234"/>
+   <Atom name="C138" type="1" charge="-0.000287"/>
+   <Atom name="C139" type="1" charge="-0.000576"/>
+   <Atom name="C140" type="1" charge="-0.015504"/>
+   <Atom name="C141" type="1" charge="-0.000775"/>
+   <Atom name="C142" type="1" charge="-0.002996"/>
+   <Atom name="C143" type="1" charge="-0.013164"/>
+   <Atom name="N001" type="2" charge="-0.018003"/>
+  </Residue>
+ </Residues>
+
+ <CoulombForce  lmax="0" coulomb14scale="1.0" 
+     mScale12="1.00" mScale13="1.00" mScale14="1.00" mScale15="1.00" mScale16="1.00"
+     pScale12="0.00" pScale13="0.00" pScale14="0.00" pScale15="0.00" pScale16="0.00" 
+     dScale12="0.00" dScale13="0.00" dScale14="0.00" dScale15="0.00" dScale16="0.00">
+   <Atom type="1"
+             c0="-0.803721"
+             />
+   <Atom type="2"
+             c0="-0.803721"
+             />
+   <Atom type="3"
+             c0="-0.803721"
+             />
+   <Atom type="4"
+             c0="-0.803721"
+             />
+   <Atom type="5"
+             c0="-0.803721"
+             />
+ </CoulombForce> 
+ <ADMPQeqForce  DampMod="3" NeutralFlag="True" SlabFlag="False" ConstQFlag="True" PbcFlag="True" 
+            >
+   <QeqAtom type="1"
+             chi="1.39198761e+02"  J="-6.51683561e-01" eta="0.08360402"  
+             />
+   <QeqAtom type="2"
+             chi="1.47411181e+02"  J="-3.42450470e-02" eta="0.08360402" 
+             />
+   <QeqAtom type="3"
+             chi="7.58342711e+01"  J="-5.62083783e-01" eta="0.08360402" 
+             />
+   <QeqAtom type="4"
+             chi="1.40533063e+02"  J="-6.33544439e-01" eta="0.08360402" 
+             />
+   <QeqAtom type="5"
+             chi="1.37525951e+02"  J="-6.63409032e-01" eta="0.08360402" 
+             />
+ </ADMPQeqForce>
+</ForceField>