Logo Search packages:      
Sourcecode: cdk version File versions  Download package

IAtomContainer org::openscience::cdk::charges::GasteigerPEPEPartialCharges::assignGasteigerPiPartialCharges ( IAtomContainer  ac,
boolean  setCharge 
) throws Exception [inline]

Main method which assigns Gasteiger partial pi charges.

Parameters:
ac AtomContainer
addCharge unused
Returns:
AtomContainer with partial charges
Exceptions:
Exception Possible Exceptions

Definition at line 110 of file GasteigerPEPEPartialCharges.java.

References org::openscience::cdk::interfaces::IAtomContainerSet::add(), org::openscience::cdk::charges::GasteigerMarsiliPartialCharges::assignGasteigerMarsiliSigmaPartialCharges(), assignPiFactors(), org::openscience::cdk::interfaces::IChemObject::clone(), fE, fS, org::openscience::cdk::tools::StructureResonanceGenerator::getAllStructures(), org::openscience::cdk::interfaces::IAtomContainer::getAtom(), org::openscience::cdk::interfaces::IAtomContainerSet::getAtomContainer(), org::openscience::cdk::interfaces::IAtomContainerSet::getAtomContainerCount(), org::openscience::cdk::interfaces::IAtomContainer::getAtomCount(), org::openscience::cdk::interfaces::IAtom::getCharge(), getElectrostaticPotentialN(), org::openscience::cdk::interfaces::IChemObject::getFlag(), org::openscience::cdk::interfaces::IAtomType::getFormalCharge(), getHyperconjugationInteractions(), getTopologicalFactors(), ISCHANGEDFC, MX_ITERATIONS, MX_RESON, org::openscience::cdk::interfaces::IAtom::setCharge(), org::openscience::cdk::interfaces::IChemObject::setFlag(), and org::openscience::cdk::charges::GasteigerMarsiliPartialCharges::setMaxGasteigerIters().

Referenced by org::openscience::cdk::qsar::descriptors::atomic::PartialPiChargeDescriptor::calculate().

                                                                                                                   {
//          logger.debug("smiles1: "+(new SmilesGenerator()).createSMILES((IMolecule) ac));
            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 gR = new StructureResonanceGenerator(true,true,true,true,false,false,MX_RESON);/*according G. should be integrated the breaking bonding*/
            IAtomContainerSet iSet = gR.getAllStructures(ac);
//          logger.debug("iset: "+iSet.getAtomContainerCount());
            
            /* detect hyperconjugation interactions */
            setHI = getHyperconjugationInteractions(ac, iSet);

            if(setHI != null) {
                  if(   setHI.getAtomContainerCount() != 0)
                        iSet.add(setHI);
//          logger.debug("isetTT: "+iSet.getAtomContainerCount());
            }
            if(iSet.getAtomContainerCount() < 2)
                  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();
//                      logger.debug("smilesClo: "+(new SmilesGenerator(ac.getBuilder())).createSMILES((IMolecule) acCloned));
                        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 k = 1 ; k < iSet.getAtomContainerCount() ; k++){
                        IAtomContainer iac = iSet.getAtomContainer(k);
//                      logger.debug("smilesITERa: "+(new SmilesGenerator(ac.getBuilder())).createSMILES((IMolecule) iac));
                        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)/* TODO- not limited*/
                                    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);
                        /*diferency 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);
                        
                        /* atom1 */
                        if(iac.getAtom(atom1).getFormalCharge() == 0){
                              if(ac.getAtom(atom1).getFormalCharge() < 0){
                                    gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = W;
                              }else{
                                    gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = -1*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] = W;
                              }else{
                                    gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = -1*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<|"+chargeT+"=c:" +charge + "+g: "+gasteigerFactors[k][STEP_SIZE * i + i + 5]);
                                    ac.getAtom(i).setCharge(chargeT);
                              }
                  }
                  
            }// iterations
            
            return ac;
            
      }


Generated by  Doxygen 1.6.0   Back to index