GasteigerPEPEPartialCharges implements IChargeCalculator { public IAtomContainer assignGasteigerPiPartialCharges(IAtomContainer ac, boolean setCharge) throws Exception { 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; 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); } StructureResonanceGenerator gR1 = new StructureResonanceGenerator(); 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(); 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); StructureResonanceGenerator gRN = new StructureResonanceGenerator(); IAtomContainerSet acSet = gRN.getContainers(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); 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++) { iSet.addAtomContainer(a.getAtomContainer(j)); } } ac = setFlags(container, ac, false); 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++) { iSet.addAtomContainer(ab.getAtomContainer(j)); } ac = setAntiFlags(container, aa, number, false); } } } 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; } 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++; } } GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges();; peoe.setMaxGasteigerIters(6); IAtomContainer acCloned; double[][] gasteigerFactors = assignPiFactors(iSet); 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); } } for (int iter = 0; iter < MX_ITERATIONS; 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) 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); 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 + ")"); double WE = Dq + fE * SumQN; logger.debug("WE : " + WE + " = Dq(" + Dq + ")+fE(" + fE + ")*SumQN(" + SumQN); int iTE = iter + 1; 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); 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; } 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); } } } logger.debug("final"); 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; } GasteigerPEPEPartialCharges(); void setMaxGasteigerIters(int iters); void setMaxResoStruc(int numbReson); int getMaxGasteigerIters(); int getMaxResoStruc(); IAtomContainer assignGasteigerPiPartialCharges(IAtomContainer ac, boolean setCharge); @Override void calculateCharges(IAtomContainer container); int getStepSize(); void setStepSize(int step); double[][] assignrPiMarsilliFactors(IAtomContainerSet setAc);  }
@Test public void testAssignGasteigerPiPartialCharges_IAtomContainer_Boolean() throws Exception { double[] testResult = {0.0, 0.0, 0.0, 0.0, 0.0}; GasteigerPEPEPartialCharges peoe = new GasteigerPEPEPartialCharges(); IAtomContainer molecule = builder.newInstance(IAtomContainer.class); molecule.addAtom(new Atom("C")); molecule.addAtom(new Atom("F")); molecule.addBond(0, 1, IBond.Order.SINGLE); addExplicitHydrogens(molecule); AtomContainerManipulator.percieveAtomTypesAndConfigureAtoms(molecule); lpcheck.saturate(molecule); peoe.assignGasteigerPiPartialCharges(molecule, true); for (int i = 0; i < molecule.getAtomCount(); i++) { Assert.assertEquals(testResult[i], molecule.getAtom(i).getCharge(), 0.01); } }