class GasteigerPEPEPartialCharges implements IChargeCalculator {

    public IAtomContainer assignGasteigerPiPartialCharges(IAtomContainer ac, boolean setCharge) throws Exception {

        // we save the aromaticity flags for the input molecule so that
        // we can add them back before we return
        boolean[] oldBondAromaticity = new boolean[ac.getBondCount()];
        boolean[] oldAtomAromaticity = new boolean[ac.getAtomCount()];
        for (int i = 0; i < ac.getAtomCount(); i++)
            oldAtomAromaticity[i] = ac.getAtom(i).getFlag(CDKConstants.ISAROMATIC);
        for (int i = 0; i < ac.getBondCount(); i++)
            oldBondAromaticity[i] = ac.getBond(i).getFlag(CDKConstants.ISAROMATIC);

        IAtomContainerSet setHI = null;

        /* 0: remove charge, and possible flag ac */
        for (int j = 0; j < ac.getAtomCount(); j++) {
            ac.getAtom(j).setCharge(0.0);
            ac.getAtom(j).setFlag(ISCHANGEDFC, false);
        }
        for (int j = 0; j < ac.getBondCount(); j++) {
            ac.getBond(j).setFlag(ISCHANGEDFC, false);
        }

        /* 1: detect resonance structure */
        StructureResonanceGenerator gR1 = new StructureResonanceGenerator();/*
                                                                             * according
                                                                             * G
                                                                             * .
                                                                             * should
                                                                             * be
                                                                             * integrated
                                                                             * the
                                                                             * breaking
                                                                             * bonding
                                                                             */
        List<IReactionProcess> reactionList1 = gR1.getReactions();
        List<IParameterReact> paramList1 = new ArrayList<IParameterReact>();
        IParameterReact param = new SetReactionCenter();
        param.setParameter(Boolean.TRUE);
        paramList1.add(param);
        HeterolyticCleavagePBReaction reactionHCPB = new HeterolyticCleavagePBReaction();
        reactionHCPB.setParameterList(paramList1);
        reactionList1.add(new SharingAnionReaction());
        Iterator<IReactionProcess> itReaction = reactionList1.iterator();
        while (itReaction.hasNext()) {
            IReactionProcess reaction = itReaction.next();
            reaction.setParameterList(paramList1);
        }
        gR1.setReactions(reactionList1);

        StructureResonanceGenerator gR2 = new StructureResonanceGenerator();/*
                                                                             * according
                                                                             * G
                                                                             * .
                                                                             * should
                                                                             * be
                                                                             * integrated
                                                                             * the
                                                                             * breaking
                                                                             * bonding
                                                                             */
        gR2.setMaximalStructures(MX_RESON);
        List<IReactionProcess> reactionList2 = gR2.getReactions();
        List<IParameterReact> paramList = new ArrayList<IParameterReact>();
        IParameterReact paramA = new SetReactionCenter();
        paramA.setParameter(Boolean.TRUE);
        paramList.add(paramA);
        reactionList2.add(new HeterolyticCleavagePBReaction());
        reactionList2.add(new SharingAnionReaction());
        itReaction = reactionList2.iterator();
        while (itReaction.hasNext()) {
            IReactionProcess reaction = itReaction.next();
            reaction.setParameterList(paramList);
        }
        gR2.setReactions(reactionList2);

        /* find resonance containers, which eliminates the repetitions */
        StructureResonanceGenerator gRN = new StructureResonanceGenerator();/*
                                                                             * according
                                                                             * G
                                                                             * .
                                                                             * should
                                                                             * be
                                                                             * integrated
                                                                             * the
                                                                             * breaking
                                                                             * bonding
                                                                             */
        IAtomContainerSet acSet = gRN.getContainers(removingFlagsAromaticity(ac));
        //		IAtomContainerSet acSet = ConjugatedPiSystemsDetector.detect(removingFlagsAromaticity(ac));

        IAtomContainerSet iSet = ac.getBuilder().newInstance(IAtomContainerSet.class);
        iSet.addAtomContainer(ac);

        if (acSet != null) for (Iterator<IAtomContainer> it = acSet.atomContainers().iterator(); it.hasNext();) {
            IAtomContainer container = it.next();
            ac = setFlags(container, ac, true);

            // Aromatic don't break its double bond homolytically
            if (Aromaticity.cdkLegacy().apply(ac))
                reactionList1.remove(reactionHCPB);
            else
                reactionList1.add(reactionHCPB);

            IAtomContainerSet a = gR1.getStructures(removingFlagsAromaticity(ac));
            if (a.getAtomContainerCount() > 1) {
                for (int j = 1; j < a.getAtomContainerCount(); j++) { // the first is already added
                    iSet.addAtomContainer(a.getAtomContainer(j));
                }
            }
            ac = setFlags(container, ac, false);

            /* processing for bonds which are not in resonance */
            for (int number = 0; number < ac.getBondCount(); number++) {
                IAtomContainer aa = setAntiFlags(container, ac, number, true);
                if (aa != null) {
                    IAtomContainerSet ab = gR2.getStructures(aa);
                    if (ab.getAtomContainerCount() > 1) for (int j = 1; j < ab.getAtomContainerCount(); j++) { // the first is already added
                            iSet.addAtomContainer(ab.getAtomContainer(j));
                        }
                    ac = setAntiFlags(container, aa, number, false);
                }
            }
        }

        /* detect hyperconjugation interactions */
        setHI = getHyperconjugationInteractions(ac, iSet);

        if (setHI != null) {
            if (setHI.getAtomContainerCount() != 0) iSet.add(setHI);
            logger.debug("setHI: " + iSet.getAtomContainerCount());
        }
        if (iSet.getAtomContainerCount() < 2) {
            for (int i = 0; i < ac.getAtomCount(); i++)
                ac.getAtom(i).setFlag(CDKConstants.ISAROMATIC, oldAtomAromaticity[i]);
            for (int i = 0; i < ac.getBondCount(); i++)
                ac.getBond(i).setFlag(CDKConstants.ISAROMATIC, oldBondAromaticity[i]);
            return ac;
        }

        /*
         * 2: search whose atoms which don't keep their formal charge and set flags
         */
        double[][] sumCharges = new double[iSet.getAtomContainerCount()][ac.getAtomCount()];
        for (int i = 1; i < iSet.getAtomContainerCount(); i++) {
            IAtomContainer iac = iSet.getAtomContainer(i);
            for (int j = 0; j < iac.getAtomCount(); j++)
                sumCharges[i][j] = iac.getAtom(j).getFormalCharge();

        }

        for (int i = 1; i < iSet.getAtomContainerCount(); i++) {
            IAtomContainer iac = iSet.getAtomContainer(i);
            int count = 0;
            for (int j = 0; j < ac.getAtomCount(); j++)
                if (count < 2) if (sumCharges[i][j] != ac.getAtom(j).getFormalCharge()) {
                    ac.getAtom(j).setFlag(ISCHANGEDFC, true);
                    iac.getAtom(j).setFlag(ISCHANGEDFC, true);
                    count++; /* TODO- error */
                }
        }

        /* 3: set sigma charge (PEOE). Initial start point */
        GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges();;
        peoe.setMaxGasteigerIters(6);
        IAtomContainer acCloned;

        double[][] gasteigerFactors = assignPiFactors(iSet);//a,b,c,deoc,chi,q

        /* 4: calculate topological weight factors Wt=fQ*fB*fA */
        double[] Wt = new double[iSet.getAtomContainerCount() - 1];
        for (int i = 1; i < iSet.getAtomContainerCount(); i++) {
            Wt[i - 1] = getTopologicalFactors(iSet.getAtomContainer(i), ac);
            logger.debug(", W:" + Wt[i - 1]);
            try {
                acCloned = (IAtomContainer) iSet.getAtomContainer(i).clone();

                acCloned = peoe.assignGasteigerMarsiliSigmaPartialCharges(acCloned, true);
                for (int j = 0; j < acCloned.getAtomCount(); j++)
                    if (iSet.getAtomContainer(i).getAtom(j).getFlag(ISCHANGEDFC)) {
                        gasteigerFactors[i][STEP_SIZE * j + j + 5] = acCloned.getAtom(j).getCharge();
                    }
            } catch (CloneNotSupportedException e) {
                throw new CDKException("Could not clone ac", e);
            }
        }

        /*
         * calculate electronegativity for changed atoms and make the difference
         * between whose atoms which change their formal charge
         */
        for (int iter = 0; iter < MX_ITERATIONS; iter++) {
            //		for (int iter = 0; iter < 1; iter++) {
            for (int k = 1; k < iSet.getAtomContainerCount(); k++) {
                IAtomContainer iac = iSet.getAtomContainer(k);
                double[] electronegativity = new double[2];
                int count = 0;
                int atom1 = 0;
                int atom2 = 0;
                for (int j = 0; j < iac.getAtomCount(); j++) {
                    if (count == 2) /*
                                     * The change of sign is product of only two
                                     * atoms, is not true
                                     */
                    break;
                    if (iac.getAtom(j).getFlag(ISCHANGEDFC)) {
                        logger.debug("Atom: " + j + ", S:" + iac.getAtom(j).getSymbol() + ", C:"
                                + iac.getAtom(j).getFormalCharge());
                        if (count == 0)
                            atom1 = j;
                        else
                            atom2 = j;

                        double q1 = gasteigerFactors[k][STEP_SIZE * j + j + 5];
                        electronegativity[count] = gasteigerFactors[k][STEP_SIZE * j + j + 2] * q1 * q1
                                + gasteigerFactors[k][STEP_SIZE * j + j + 1] * q1
                                + gasteigerFactors[k][STEP_SIZE * j + j];
                        logger.debug("e:" + electronegativity[count] + ",q1: " + q1 + ", c:"
                                + gasteigerFactors[k][STEP_SIZE * j + j + 2] + ", b:"
                                + gasteigerFactors[k][STEP_SIZE * j + j + 1] + ", a:"
                                + gasteigerFactors[k][STEP_SIZE * j + j]);
                        count++;
                    }
                }
                logger.debug("Atom1:" + atom1 + ",Atom2:" + atom2);
                /* difference of electronegativity 1 lower */
                double max1 = Math.max(electronegativity[0], electronegativity[1]);
                double min1 = Math.min(electronegativity[0], electronegativity[1]);
                double DX = 1.0;
                if (electronegativity[0] < electronegativity[1])
                    DX = gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 3];
                else
                    DX = gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 3];

                double Dq = (max1 - min1) / DX;
                logger.debug("Dq : " + Dq + " = (" + max1 + "-" + min1 + ")/" + DX);
                double epN1 = getElectrostaticPotentialN(iac, atom1, gasteigerFactors[k]);
                double epN2 = getElectrostaticPotentialN(iac, atom2, gasteigerFactors[k]);
                double SumQN = Math.abs(epN1 - epN2);
                logger.debug("sum(" + SumQN + ") = (" + epN1 + ") - (" + epN2 + ")");

                /* electronic weight */
                double WE = Dq + fE * SumQN;
                logger.debug("WE : " + WE + " = Dq(" + Dq + ")+fE(" + fE + ")*SumQN(" + SumQN);
                int iTE = iter + 1;

                /* total topological */
                double W = WE * Wt[k - 1] * fS / (iTE);
                logger.debug("W : " + W + " = WE(" + WE + ")*Wt(" + Wt[k - 1] + ")*fS(" + fS + ")/iter(" + iTE
                        + "), atoms: " + atom1 + ", " + atom2);

                /* iac == new structure, ac == old structure */
                /* atom1 */
                if (iac.getAtom(atom1).getFormalCharge() == 0) {
                    if (ac.getAtom(atom1).getFormalCharge() < 0) {
                        gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = -1 * W;
                    } else {
                        gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = W;
                    }
                } else if (iac.getAtom(atom1).getFormalCharge() > 0) {
                    gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = W;
                } else {
                    gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = -1 * W;
                }
                /* atom2 */
                if (iac.getAtom(atom2).getFormalCharge() == 0) {
                    if (ac.getAtom(atom2).getFormalCharge() < 0) {
                        gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = -1 * W;
                    } else {
                        gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = W;
                    }
                } else if (iac.getAtom(atom2).getFormalCharge() > 0) {
                    gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = W;
                } else {
                    gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = -1 * W;
                }

            }
            for (int k = 1; k < iSet.getAtomContainerCount(); k++) {

                for (int i = 0; i < ac.getAtomCount(); i++)
                    if (iSet.getAtomContainer(k).getAtom(i).getFlag(ISCHANGEDFC)) {
                        double charge = ac.getAtom(i).getCharge();
                        double chargeT = 0.0;
                        chargeT = charge + gasteigerFactors[k][STEP_SIZE * i + i + 5];
                        logger.debug("i<|" + ac.getAtom(i).getSymbol() + ", " + chargeT + "=c:" + charge + "+g: "
                                + gasteigerFactors[k][STEP_SIZE * i + i + 5]);
                        ac.getAtom(i).setCharge(chargeT);
                    }
            }

        }// iterations
        logger.debug("final");

        // before getting back we should set back the aromatic flags
        for (int i = 0; i < ac.getAtomCount(); i++)
            ac.getAtom(i).setFlag(CDKConstants.ISAROMATIC, oldAtomAromaticity[i]);
        for (int i = 0; i < ac.getBondCount(); i++)
            ac.getBond(i).setFlag(CDKConstants.ISAROMATIC, oldBondAromaticity[i]);

        return ac;

    }

    public  GasteigerPEPEPartialCharges();

}