← Back to branch summary

~ubuntu-branches/ubuntu/maverick/cdk/maverick

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/charges/GasteigerPEPEPartialCharges.java

  • Committer: Bazaar Package Importer
  • Author(s): Paul Cager
  • Date: 2008-04-09 21:17:53 UTC
  • Revision ID: james.westby@ubuntu.com-20080409211753-46lmjw5z8mx5pd8d
Tags: upstream-1.0.2
Import upstream version 1.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/org/openscience/cdk/charges/GasteigerPEPEPartialCharges.java
 
1
/*  $RCSfile$
 
2
 *  $Author: miguelrojasch $
 
3
 *  $Date: 2006-05-11 10:17:36 +0200 (Do, 11 Mai 2006) $
 
4
 *  $Revision: 6217 $
 
5
 *
 
6
 *  Copyright (C) 2004-2007  Miguel Rojas <miguel.rojas@uni-koeln.de>
 
7
 *
 
8
 *  Contact: cdk-devel@list.sourceforge.net
 
9
 *
 
10
 *  This program is free software; you can redistribute it and/or
 
11
 *  modify it under the terms of the GNU Lesser General Public License
 
12
 *  as published by the Free Software Foundation; either version 2.1
 
13
 *  of the License, or (at your option) any later version.
 
14
 *
 
15
 *  This program is distributed in the hope that it will be useful,
 
16
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
 *  GNU Lesser General Public License for more details.
 
19
 *
 
20
 *  You should have received a copy of the GNU Lesser General Public License
 
21
 *  along with this program; if not, write to the Free Software
 
22
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 
23
 */
 
24
package org.openscience.cdk.charges;
 
25
 
 
26
import org.openscience.cdk.CDKConstants;
 
27
import org.openscience.cdk.aromaticity.HueckelAromaticityDetector;
 
28
import org.openscience.cdk.config.AtomTypeFactory;
 
29
import org.openscience.cdk.exception.CDKException;
 
30
import org.openscience.cdk.interfaces.*;
 
31
import org.openscience.cdk.reaction.IReactionProcess;
 
32
import org.openscience.cdk.reaction.type.BreakingBondReaction;
 
33
import org.openscience.cdk.reaction.type.HyperconjugationReaction;
 
34
import org.openscience.cdk.tools.LoggingTool;
 
35
import org.openscience.cdk.tools.StructureResonanceGenerator;
 
36
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
 
37
 
 
38
import java.io.IOException;
 
39
 
 
40
/**
 
41
 * <p>The calculation of the Gasteiger (PEPE) partial charges is based on 
 
42
 * {@cdk.cite Saller85}. This class doesn't implement the original method of the Marsili but the 
 
43
 * method based on H. Saller which is described from Petra manual version 2.6</p> 
 
44
 * <p>They are calculated by generating all valence bond(resonance) structures
 
45
 * for this system and then weighting them on the basis of pi-orbital electronegativies
 
46
 * and formal considerations based on PEPE (Partial Equalization of pi-electronegativity).</p>
 
47
 * 
 
48
 * @author      Miguel Rojas
 
49
 * 
 
50
 * @cdk.module  charges
 
51
 * @cdk.created 2006-05-14
 
52
 * @cdk.keyword partial atomic charges
 
53
 * @cdk.keyword charge distribution
 
54
 * @cdk.keyword electronegativities, partial equalization of orbital
 
55
 * @cdk.keyword PEPE
 
56
 * @see GasteigerMarsiliPartialCharges
 
57
 */
 
58
public class GasteigerPEPEPartialCharges {
 
59
        /** max iterations */
 
60
        private double MX_ITERATIONS = 8;
 
61
        /** max number of resonance structures to be searched*/
 
62
        private int MX_RESON = 50;
 
63
        private int STEP_SIZE = 5;
 
64
        private AtomTypeFactory factory;
 
65
        /** Flag is set if the formal charge of a chemobject is changed due to resonance.*/
 
66
        private static int ISCHANGEDFC = 0;
 
67
        
 
68
        /** Corresponds an empirical influence between the electrostatic potential and
 
69
         * the neighbours.*/
 
70
        private double fE = 1.1;/*1.1*/
 
71
        /** Scalle factor which makes same heavay for all structures*/
 
72
        private double fS = 0.37;
 
73
        
 
74
        
 
75
        
 
76
        private LoggingTool logger = new LoggingTool(GasteigerPEPEPartialCharges.class);
 
77
 
 
78
        
 
79
        /**
 
80
         *  Constructor for the GasteigerPEPEPartialCharges object
 
81
         */
 
82
        public GasteigerPEPEPartialCharges() { }
 
83
        /**
 
84
         *  Sets the maxGasteigerIters attribute of the GasteigerPEPEPartialCharges
 
85
         *  object
 
86
         *
 
87
         *@param  iters  The new maxGasteigerIters value
 
88
         */
 
89
        public void setMaxGasteigerIters(double iters) {
 
90
                MX_ITERATIONS = iters;
 
91
        }
 
92
        /**
 
93
         *  Sets the maximum resonance structures to be searched
 
94
         *
 
95
         *@param  numbReson  The number of resonance Structures to be searched
 
96
         */
 
97
        public void setMaxResoStruc(int numbReson) {
 
98
                MX_RESON = numbReson;
 
99
        }
 
100
        
 
101
        /**
 
102
         *  Main method which assigns Gasteiger partial pi charges. 
 
103
         *  
 
104
         *
 
105
         *@param  ac             AtomContainer
 
106
         *@param  addCharge      unused
 
107
         *@return                AtomContainer with partial charges
 
108
         *@exception  Exception  Possible Exceptions
 
109
         */
 
110
        public IAtomContainer assignGasteigerPiPartialCharges(IAtomContainer ac, boolean setCharge) throws Exception {
 
111
//              logger.debug("smiles1: "+(new SmilesGenerator()).createSMILES((IMolecule) ac));
 
112
                IAtomContainerSet setHI = null;
 
113
                
 
114
                /*0: remove charge, and possible flag ac*/
 
115
                for(int j = 0 ; j < ac.getAtomCount(); j++){
 
116
                        ac.getAtom(j).setCharge(0.0);
 
117
                        ac.getAtom(j).setFlag(ISCHANGEDFC, false);
 
118
                }
 
119
                for(int j = 0 ; j < ac.getBondCount(); j++){
 
120
                        ac.getBond(j).setFlag(ISCHANGEDFC, false);
 
121
                }
 
122
                
 
123
                /*1: detect resonance structure*/
 
124
                StructureResonanceGenerator gR = new StructureResonanceGenerator(true,true,true,true,false,false,MX_RESON);/*according G. should be integrated the breaking bonding*/
 
125
                IAtomContainerSet iSet = gR.getAllStructures(ac);
 
126
//              logger.debug("iset: "+iSet.getAtomContainerCount());
 
127
                
 
128
                /* detect hyperconjugation interactions */
 
129
                setHI = getHyperconjugationInteractions(ac, iSet);
 
130
 
 
131
                if(setHI != null) {
 
132
                        if(     setHI.getAtomContainerCount() != 0)
 
133
                                iSet.add(setHI);
 
134
//              logger.debug("isetTT: "+iSet.getAtomContainerCount());
 
135
                }
 
136
                if(iSet.getAtomContainerCount() < 2)
 
137
                        return ac;
 
138
                
 
139
                /*2: search whose atoms which don't keep their formal charge and set flags*/
 
140
                double[][] sumCharges = new double[iSet.getAtomContainerCount()][ac.getAtomCount( )];
 
141
                for(int i = 1; i < iSet.getAtomContainerCount() ; i++){
 
142
                        IAtomContainer iac = iSet.getAtomContainer(i);
 
143
                        for(int j = 0 ; j < iac.getAtomCount(); j++)
 
144
                                sumCharges[i][j] = iac.getAtom(j).getFormalCharge();
 
145
                        
 
146
                }
 
147
                for(int i = 1; i < iSet.getAtomContainerCount() ; i++){
 
148
                        IAtomContainer iac = iSet.getAtomContainer(i);
 
149
                        int count = 0;
 
150
                        for(int j = 0 ; j < ac.getAtomCount(); j++)
 
151
                                if(count < 2)
 
152
                                if(sumCharges[i][j] != ac.getAtom(j).getFormalCharge()){
 
153
                                        ac.getAtom(j).setFlag(ISCHANGEDFC, true);
 
154
                                        iac.getAtom(j).setFlag(ISCHANGEDFC, true);
 
155
                                        count++; /* TODO- error*/
 
156
                                }
 
157
                }
 
158
 
 
159
                /*3: set sigma charge (PEOE). Initial start point*/
 
160
                GasteigerMarsiliPartialCharges peoe = new GasteigerMarsiliPartialCharges();;
 
161
                peoe.setMaxGasteigerIters(6);
 
162
                IAtomContainer acCloned;
 
163
                
 
164
 
 
165
                double[][] gasteigerFactors = assignPiFactors(iSet);//a,b,c,deoc,chi,q
 
166
                
 
167
                /*4: calculate topological weight factors Wt=fQ*fB*fA*/
 
168
                double[] Wt = new double[iSet.getAtomContainerCount()-1];
 
169
                for(int i = 1; i < iSet.getAtomContainerCount() ; i++){
 
170
                        Wt[i-1]= getTopologicalFactors(iSet.getAtomContainer(i),ac);
 
171
//                      logger.debug(", W:"+Wt[i-1]);
 
172
                        try {
 
173
                                acCloned = (IAtomContainer)iSet.getAtomContainer(i).clone();
 
174
//                              logger.debug("smilesClo: "+(new SmilesGenerator(ac.getBuilder())).createSMILES((IMolecule) acCloned));
 
175
                                acCloned = peoe.assignGasteigerMarsiliSigmaPartialCharges(acCloned, true);
 
176
                                for(int j = 0 ; j<acCloned.getAtomCount(); j++)
 
177
                                        if(iSet.getAtomContainer(i).getAtom(j).getFlag(ISCHANGEDFC)){
 
178
                                                gasteigerFactors[i][STEP_SIZE * j + j + 5] = acCloned.getAtom(j).getCharge(); 
 
179
                                        }
 
180
                        } catch (CloneNotSupportedException e) {
 
181
                                throw new CDKException("Could not clone ac", e);
 
182
                        }
 
183
                }
 
184
                
 
185
                /*calculate electronegativity for changed atoms and make the difference between whose
 
186
                 * atoms which change their formal charge*/
 
187
                for (int iter = 0; iter < MX_ITERATIONS; iter++) {
 
188
                        for(int k = 1 ; k < iSet.getAtomContainerCount() ; k++){
 
189
                                IAtomContainer iac = iSet.getAtomContainer(k);
 
190
//                              logger.debug("smilesITERa: "+(new SmilesGenerator(ac.getBuilder())).createSMILES((IMolecule) iac));
 
191
                                double[] electronegativity = new double[2];
 
192
                                int count = 0;
 
193
                                int atom1 = 0;
 
194
                                int atom2 = 0;
 
195
                                for (int j = 0; j < iac.getAtomCount(); j++) {
 
196
                                        if(count == 2)/* TODO- not limited*/
 
197
                                                break;
 
198
                                        if(iac.getAtom(j).getFlag(ISCHANGEDFC)){
 
199
//                                              logger.debug("Atom: "+j+", S:"+iac.getAtom(j).getSymbol()+", C:"+iac.getAtom(j).getFormalCharge());
 
200
                                                if(count == 0)
 
201
                                                        atom1 = j;
 
202
                                                else 
 
203
                                                        atom2 = j;
 
204
                                                
 
205
                                                double q1 = gasteigerFactors[k][STEP_SIZE * j + j + 5];
 
206
                                                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];
 
207
//                                              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]);
 
208
                                                count++;
 
209
                                        }
 
210
                                        
 
211
                                }
 
212
//                              logger.debug("Atom1:"+atom1+",Atom2:"+atom2);
 
213
                                /*diferency of electronegativity 1 lower*/
 
214
                                double max1 = Math.max(electronegativity[0], electronegativity[1]);
 
215
                                double min1 = Math.min(electronegativity[0], electronegativity[1]);
 
216
                                double DX = 1.0;
 
217
                                if(electronegativity[0] < electronegativity[1])
 
218
                                        DX = gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 3];
 
219
                                else
 
220
                                        DX = gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 3];
 
221
                                        
 
222
                                double Dq = (max1-min1)/DX;
 
223
//                                      logger.debug("Dq : "+Dq+ " = ("+ max1+"-"+min1+")/"+DX);
 
224
                                double epN1 = getElectrostaticPotentialN(iac,atom1,gasteigerFactors[k]);
 
225
                                double epN2 = getElectrostaticPotentialN(iac,atom2,gasteigerFactors[k]);
 
226
                                double SumQN = Math.abs(epN1 - epN2);
 
227
//                                      logger.debug("sum("+SumQN+") = ("+epN1+") - ("+epN2+")");
 
228
                                
 
229
                                /* electronic weight*/
 
230
                                double WE = Dq + fE*SumQN;
 
231
//                                      logger.debug("WE : "+WE+" = Dq("+Dq+")+fE("+fE+")*SumQN("+SumQN);
 
232
                                int iTE = iter+1;
 
233
                                
 
234
                                /* total topological*/
 
235
                                double W = WE*Wt[k-1]*fS/(iTE);
 
236
//                                      logger.debug("W : "+W+" = WE("+WE+")*Wt("+Wt[k-1]+")*fS("+fS+")/iter("+iTE+"), atoms: "+atom1+", "+atom2);
 
237
                                
 
238
                                /* atom1 */
 
239
                                if(iac.getAtom(atom1).getFormalCharge() == 0){
 
240
                                        if(ac.getAtom(atom1).getFormalCharge() < 0){
 
241
                                                gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = W;
 
242
                                        }else{
 
243
                                                gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = -1*W;
 
244
                                        }
 
245
                                }else if(iac.getAtom(atom1).getFormalCharge() > 0){
 
246
                                        gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = W;
 
247
                                }else{
 
248
                                        gasteigerFactors[k][STEP_SIZE * atom1 + atom1 + 5] = -1*W;
 
249
                                }
 
250
                                /* atom2*/
 
251
                                if(iac.getAtom(atom2).getFormalCharge() == 0){
 
252
                                        if(ac.getAtom(atom2).getFormalCharge() < 0){
 
253
                                                gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = W;
 
254
                                        }else{
 
255
                                                gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = -1*W;
 
256
                                        }
 
257
                                }else if(iac.getAtom(atom2).getFormalCharge() > 0){
 
258
                                        gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = W;
 
259
                                }else{
 
260
                                        gasteigerFactors[k][STEP_SIZE * atom2 + atom2 + 5] = -1*W;
 
261
                                }
 
262
                                
 
263
                        }
 
264
                        for(int k = 1 ; k < iSet.getAtomContainerCount() ; k++){
 
265
                                for (int i = 0; i < ac.getAtomCount(); i++) 
 
266
                                        if(iSet.getAtomContainer(k).getAtom(i).getFlag(ISCHANGEDFC)){
 
267
                                                double charge = ac.getAtom(i).getCharge();
 
268
                                                double chargeT = 0.0;
 
269
                                                chargeT = charge + gasteigerFactors[k][STEP_SIZE * i + i + 5];
 
270
//                                              logger.debug("i<|"+chargeT+"=c:" +charge + "+g: "+gasteigerFactors[k][STEP_SIZE * i + i + 5]);
 
271
                                                ac.getAtom(i).setCharge(chargeT);
 
272
                                        }
 
273
                        }
 
274
                        
 
275
                }// iterations
 
276
                
 
277
                return ac;
 
278
                
 
279
        }
 
280
        /**
 
281
         * get the possibles structures after a hyperconjugation interactions for bonds wich
 
282
         * not belong any resonance structure.
 
283
         * 
 
284
         * @param ac IAtomContainer
 
285
         * @return IAtomContainerSet
 
286
         * @throws CDKException 
 
287
         * @throws ClassNotFoundException 
 
288
         * @throws IOException 
 
289
         */
 
290
        private IAtomContainerSet getHyperconjugationInteractions(IAtomContainer ac, IAtomContainerSet iSet) throws IOException, ClassNotFoundException, CDKException {
 
291
                IAtomContainerSet set = ac.getBuilder().newAtomContainerSet();
 
292
        IReactionProcess type = new BreakingBondReaction();
 
293
        cleanFlagReactiveCenter(ac);
 
294
        boolean found = false; /* control obtained containers */
 
295
                IMoleculeSet setOfReactants = ac.getBuilder().newMoleculeSet();
 
296
                /* search of reactive center.*/
 
297
                out:
 
298
                for(int i = 0 ; i < ac.getBondCount() ; i++){
 
299
                        if(ac.getBond(i).getOrder() > 1 ){
 
300
                                for(int j = 0 ; j < iSet.getAtomContainerCount(); j++){
 
301
                                IAtomContainer ati = iSet.getAtomContainer(j);
 
302
                                if(!ati.equals(ac))
 
303
                                for(int k = 0; k < ati.getBondCount(); k++){
 
304
                                        IAtom a0 = ati.getBond(k).getAtom(0);
 
305
                                        IAtom a1 = ati.getBond(k).getAtom(1);
 
306
                                        if(!a0.getSymbol().equals("H") || !a1.getSymbol().equals("H"))
 
307
                                        if((a0.getID().equals(ac.getBond(i).getAtom(0).getID()) &&
 
308
                                                        a1.getID().equals(ac.getBond(i).getAtom(1).getID())) ||
 
309
                                                        (a1.getID().equals(ac.getBond(i).getAtom(0).getID()) &&
 
310
                                                                a0.getID().equals(ac.getBond(i).getAtom(1).getID()))){
 
311
                                                if(a0.getFormalCharge() != 0 || a1.getFormalCharge() != 0)
 
312
                                                        continue out;
 
313
                                        }
 
314
                                }
 
315
                                }
 
316
                                ac.getBond(i).getAtom(0).setFlag(CDKConstants.REACTIVE_CENTER,true);
 
317
                                ac.getBond(i).getAtom(1).setFlag(CDKConstants.REACTIVE_CENTER,true);
 
318
                                ac.getBond(i).setFlag(CDKConstants.REACTIVE_CENTER,true);
 
319
                                found = true;
 
320
                        }
 
321
                }
 
322
                if(!found)
 
323
                        return null;
 
324
                
 
325
                
 
326
                setOfReactants.addMolecule((IMolecule) ac);
 
327
                Object[] params = {Boolean.TRUE};
 
328
                
 
329
                type.setParameters(params);
 
330
                IReactionSet setOfReactions = type.initiate(setOfReactants, null);
 
331
        for(int i = 0; i < setOfReactions.getReactionCount(); i++){
 
332
                type = new HyperconjugationReaction();
 
333
                IMoleculeSet setOfM2 = ac.getBuilder().newMoleculeSet();
 
334
                IMolecule mol= setOfReactions.getReaction(i).getProducts().getMolecule(0);
 
335
                for(int k = 0; k < mol.getBondCount(); k++){
 
336
                        mol.getBond(k).setFlag(CDKConstants.REACTIVE_CENTER,false);
 
337
                        mol.getBond(k).getAtom(0).setFlag(CDKConstants.REACTIVE_CENTER,false);
 
338
                        mol.getBond(k).getAtom(1).setFlag(CDKConstants.REACTIVE_CENTER,false);
 
339
                }
 
340
                setOfM2.addMolecule((IMolecule) mol);
 
341
                Object[] params2 = {Boolean.FALSE};
 
342
                        type.setParameters(params2);
 
343
                        IReactionSet setOfReactions2 = type.initiate(setOfM2, null);
 
344
                        if(setOfReactions2.getReactionCount() > 0){
 
345
                                
 
346
                        IMolecule react = setOfReactions2.getReaction(0).getReactants().getMolecule(0);
 
347
 
 
348
                        set.addAtomContainer(react);
 
349
                        }
 
350
        }
 
351
 
 
352
                return set;
 
353
        }
 
354
        /**
 
355
         * get the electrostatic potential of the neighbours of a atom.
 
356
         *  
 
357
         * @param ac   The IAtomContainer to study
 
358
         * @param ds 
 
359
         * @param atom1 The position of the IAtom to study
 
360
         * @return     The sum of electrostatic potential of the neighbours
 
361
         */
 
362
        private double getElectrostaticPotentialN(IAtomContainer ac, int atom1, double[] ds) {
 
363
                
 
364
//              double CoulombForceConstant = 1/(4*Math.PI*8.81/*Math.pow(10, -12)*/);
 
365
                double CoulombForceConstant = 0.048;
 
366
                double sum = 0.0;
 
367
                try {
 
368
                        if (factory == null) 
 
369
                factory = AtomTypeFactory.getInstance(
 
370
                    "org/openscience/cdk/config/data/jmol_atomtypes.txt", 
 
371
                    ac.getBuilder()
 
372
                );
 
373
 
 
374
                        java.util.List atoms = ac.getConnectedAtomsList(ac.getAtom(atom1));
 
375
                        for(int i = 0 ; i < atoms.size() ; i++){
 
376
                                double covalentradius = 0;
 
377
                    String symbol = ((IAtom)atoms.get(i)).getSymbol();
 
378
                    IAtomType type = factory.getAtomType(symbol);
 
379
                    covalentradius = type.getCovalentRadius();
 
380
 
 
381
                    double charge = ds[STEP_SIZE * atom1 + atom1 + 5];
 
382
                                double sumI = CoulombForceConstant*charge/(covalentradius*covalentradius);
 
383
//                              logger.debug("sum_("+sumI+") = CFC("+CoulombForceConstant+")*charge("+charge+"/ret("+covalentradius);
 
384
                                sum += sumI;
 
385
                        }
 
386
                } catch (CDKException e) {
 
387
            logger.debug(e);
 
388
        }
 
389
                
 
390
                return sum;
 
391
        }
 
392
 
 
393
 
 
394
        /**
 
395
         * get the topological weight factor for each atomContainer
 
396
         * 
 
397
         * @param atomContainer  The IAtomContainer to study.
 
398
         * @param ac             The IAtomContainer to study.
 
399
         * @return The value
 
400
         */
 
401
        private double getTopologicalFactors(IAtomContainer atomContainer,IAtomContainer ac) {
 
402
                /*factor for separation of charge*/
 
403
                int totalNCharge1 = AtomContainerManipulator.getTotalNegativeFormalCharge(atomContainer);
 
404
                int totalPCharge1 = AtomContainerManipulator.getTotalPositiveFormalCharge(atomContainer);
 
405
                
 
406
                double fQ = 1.0;
 
407
                if(totalNCharge1 != 0.0){
 
408
                        fQ = 0.5;
 
409
                        for(int i = 0; i < atomContainer.getBondCount(); i++){
 
410
                                IBond bond = atomContainer.getBond(i);
 
411
                                if(bond.getAtom(0).getFormalCharge() != 0.0 && bond.getAtom(1).getFormalCharge() != 0.0){
 
412
                                        fQ = 0.25;
 
413
                                        break;
 
414
                                }
 
415
                        }
 
416
                }
 
417
                /*factor, if the number of covalents bonds is decreased*/
 
418
                double fB = 1.0;
 
419
                
 
420
                int numBond1 = 0;
 
421
                int numBond2 = 0;
 
422
        for (int i = 0; i < atomContainer.getBondCount(); i++) {
 
423
            if (atomContainer.getBond(i).getOrder() == 2.0)
 
424
                numBond1 += 1;
 
425
            if (ac.getBond(i).getOrder() == 2.0)
 
426
                numBond2 += 1;
 
427
        }
 
428
        
 
429
        if(numBond1 </*>*/ numBond2)
 
430
                        fB = 0.8;
 
431
        
 
432
                double fPlus = 1.0;
 
433
        if(totalNCharge1 == 0.0 && totalPCharge1 == 0.0 )
 
434
                fPlus = 0.1;
 
435
        
 
436
        
 
437
        /*aromatic*/
 
438
        double fA = 1.0;
 
439
        try {
 
440
                        if(HueckelAromaticityDetector.detectAromaticity(ac))
 
441
                                if(!HueckelAromaticityDetector.detectAromaticity(atomContainer))
 
442
                                                fA = 0.3;
 
443
                } catch (CDKException e) {
 
444
                        e.printStackTrace();
 
445
                }
 
446
//              logger.debug("return "+fQ*fB*fPlus*fA+"= sp:"+fQ+", dc:"+fB+", fPlus:"+fPlus+", fA:"+fA);
 
447
                
 
448
                return fQ*fB*fPlus*fA;
 
449
        }
 
450
 
 
451
 
 
452
        /**
 
453
         *  Get the StepSize attribute of the GasteigerMarsiliPartialCharges
 
454
         *  object
 
455
         *
 
456
         *@return STEP_SIZE
 
457
         */
 
458
        public int getStepSize(){
 
459
                return STEP_SIZE;
 
460
        }
 
461
        
 
462
 
 
463
        /**
 
464
         * Method which stores and assigns the factors a,b,c and CHI+
 
465
         *
 
466
         * @return     Array of doubles [a1,b1,c1,denom1,chi1,q1...an,bn,cn...] 1:Atom 1-n in AtomContainer
 
467
         */
 
468
        private double[][] assignPiFactors(IAtomContainerSet setAc) {
 
469
                //a,b,c,denom,chi,q
 
470
                double[][] gasteigerFactors = new double[setAc.getAtomContainerCount()][(setAc.getAtomContainer(0).getAtomCount() * (STEP_SIZE+1))];
 
471
                String AtomSymbol = "";
 
472
                double[] factors = new double[]{0.0, 0.0, 0.0};
 
473
                for( int k = 1 ; k < setAc.getAtomContainerCount(); k ++){
 
474
                        IAtomContainer ac = setAc.getAtomContainer(k);
 
475
                        for (int i = 0; i < ac.getAtomCount(); i++) {
 
476
                                factors[0] = 0.0;
 
477
                                factors[1] = 0.0;
 
478
                                factors[2] = 0.0;
 
479
                                AtomSymbol = ac.getAtom(i).getSymbol();
 
480
                                if (AtomSymbol.equals("H")) {
 
481
                                        factors[0] = 0.0;
 
482
                                        factors[1] = 0.0;
 
483
                                        factors[2] = 0.0;
 
484
                                } else if (AtomSymbol.equals("C")) {/*
 
485
                                        if(ac.getAtom(i).getFlag(ISCHANGEDFC))*/{
 
486
                                                factors[0] = 5.60;
 
487
                                                factors[1] = 8.93;
 
488
                                                factors[2] = 2.94;
 
489
                                        }
 
490
                                } else if (AtomSymbol.equals("O")) {
 
491
                                                if(ac.getMaximumBondOrder(ac.getAtom(i)) == 1){
 
492
                                                        factors[0] = 10.0;
 
493
                                                        factors[1] = 13.86;
 
494
                                                        factors[2] = 9.68;
 
495
                                                }else {
 
496
                                                        factors[0] = 7.91;
 
497
                                                        factors[1] = 14.76;
 
498
                                                        factors[2] = 6.85;
 
499
                                                }
 
500
                                } else if (AtomSymbol.equals("N")) {
 
501
                                        if(ac.getMaximumBondOrder(ac.getAtom(i)) != 1){
 
502
                                                factors[0] = 7.95;/*7.95*/
 
503
                                                factors[1] = 9.73;/*9.73*/
 
504
                                                factors[2] = 2.67;/*2.67*/
 
505
                                        }else {
 
506
                                                factors[0] = 4.54;/*4.54*//*5.5*/
 
507
                                                factors[1] = 11.86;/*11.86*//*10.86*/
 
508
                                                factors[2] = 7.32;/*7.32*//*7.99*/
 
509
                                        }
 
510
                                } else if (AtomSymbol.equals("S")) {
 
511
                                        if(ac.getMaximumBondOrder(ac.getAtom(i)) == 1){
 
512
                                                factors[0] = 7.73;
 
513
                                                factors[1] = 8.16;
 
514
                                                factors[2] = 1.81;
 
515
                                        }else {
 
516
                                                factors[0] = 6.60;
 
517
                                                factors[1] = 10.32;
 
518
                                                factors[2] = 3.72;
 
519
                                        }
 
520
                                } else if (AtomSymbol.equals("F")) {
 
521
                                        factors[0] = 7.34;
 
522
                                        factors[1] = 13.86;
 
523
                                        factors[2] = 9.68;
 
524
                                } else if (AtomSymbol.equals("Cl")) {
 
525
                                        factors[0] = 6.50;
 
526
                                        factors[1] = 11.02;
 
527
                                        factors[2] = 4.52;
 
528
                                } else if (AtomSymbol.equals("Br")) {
 
529
                                        factors[0] = 5.20;
 
530
                                        factors[1] = 9.68;
 
531
                                        factors[2] = 4.48;
 
532
                                } else if (AtomSymbol.equals("I")) {
 
533
                                        factors[0] = 4.95;
 
534
                                        factors[1] = 8.81;
 
535
                                        factors[2] = 3.86;
 
536
                                }
 
537
                        
 
538
                                gasteigerFactors[k][STEP_SIZE * i + i] = factors[0];
 
539
                                gasteigerFactors[k][STEP_SIZE * i + i + 1] = factors[1];
 
540
                                gasteigerFactors[k][STEP_SIZE * i + i + 2] = factors[2];
 
541
                                gasteigerFactors[k][STEP_SIZE * i + i + 5] = ac.getAtom(i).getCharge();
 
542
        
 
543
                                if (factors[0] == 0 && factors[1] == 0 && factors[2] == 0) {
 
544
                                        gasteigerFactors[k][STEP_SIZE * i + i + 3] = 1;
 
545
                                } else {
 
546
                                        gasteigerFactors[k][STEP_SIZE * i + i + 3] = factors[0] + factors[1] + factors[2];
 
547
                                }
 
548
                        }
 
549
                }
 
550
                
 
551
 
 
552
                return gasteigerFactors;
 
553
        }
 
554
        /**
 
555
         *  Method which stores and assigns the factors a,b,c and CHI+
 
556
         *
 
557
         *@return     Array of doubles [a1,b1,c1,denom1,chi1,q1...an,bn,cn...] 1:Atom 1-n in AtomContainer
 
558
         */
 
559
        public double[][] assignrPiMarsilliFactors(IAtomContainerSet setAc) {
 
560
                //a,b,c,denom,chi,q
 
561
                double[][] gasteigerFactors = new double[setAc.getAtomContainerCount()][(setAc.getAtomContainer(0).getAtomCount() * (STEP_SIZE+1))];
 
562
                String AtomSymbol = "";
 
563
                double[] factors = new double[]{0.0, 0.0, 0.0};
 
564
                for( int k = 1 ; k < setAc.getAtomContainerCount(); k ++){
 
565
                        IAtomContainer ac = setAc.getAtomContainer(k);
 
566
                        
 
567
                        for (int i = 0; i < ac.getAtomCount(); i++) {
 
568
                                factors[0] = 0.0;
 
569
                                factors[1] = 0.0;
 
570
                                factors[2] = 0.0;
 
571
                                AtomSymbol = ac.getAtom(i).getSymbol();
 
572
                                if (AtomSymbol.equals("H")) {
 
573
                                        factors[0] = 0.0;
 
574
                                        factors[1] = 0.0;
 
575
                                        factors[2] = 0.0;
 
576
                                } else if (AtomSymbol.equals("C")) {
 
577
                                                factors[0] = 5.98;/*5.98-5.60*/
 
578
                                                factors[1] = 7.93;/*7.93-8.93*/
 
579
                                                factors[2] = 1.94;
 
580
                                } else if (AtomSymbol.equals("O")) {
 
581
                                                if(ac.getMaximumBondOrder(ac.getAtom(i)) > 1){
 
582
                                                        factors[0] = 11.2;/*11.2-10.0*/
 
583
                                                        factors[1] = 13.24;/*13.24-13.86*/
 
584
                                                        factors[2] = 9.68;
 
585
                                                }else {
 
586
                                                        factors[0] = 7.91;
 
587
                                                        factors[1] = 14.76;
 
588
                                                        factors[2] = 6.85;
 
589
                                                }
 
590
                                } else if (AtomSymbol.equals("N")) {
 
591
                                        if(ac.getMaximumBondOrder(ac.getAtom(i))  > 1){
 
592
                                                
 
593
                                                factors[0] = 8.95;/*7.95*/
 
594
                                                factors[1] = 9.73;/*9.73*/
 
595
                                                factors[2] = 2.67;/*2.67*/
 
596
                                        }else {
 
597
                                                factors[0] = 4.54;
 
598
                                                factors[1] = 11.86;
 
599
                                                factors[2] = 7.32;
 
600
                                        }
 
601
                                } else if (AtomSymbol.equals("P")) {// <--No correct
 
602
                                        if(ac.getMaximumBondOrder(ac.getAtom(i))  > 1){
 
603
                                                factors[0] = 10.73;// <--No correct
 
604
                                                factors[1] = 11.16;// <--No correct
 
605
                                                factors[2] = 6.81;// <--No correct
 
606
                                        }else {
 
607
                                                factors[0] = 9.60;// <--No correct
 
608
                                                factors[1] = 13.32;// <--No correct
 
609
                                                factors[2] = 2.72;// <--No correct
 
610
                                        }
 
611
                                } else if (AtomSymbol.equals("S")) {
 
612
                                        if(ac.getMaximumBondOrder(ac.getAtom(i))  > 1){
 
613
                                                
 
614
                                                factors[0] = 7.73;
 
615
                                                factors[1] = 8.16;
 
616
                                                factors[2] = 1.81;
 
617
                                        }else {
 
618
                                                factors[0] = 6.60;
 
619
                                                factors[1] = 10.32;
 
620
                                                factors[2] = 3.72;
 
621
                                        }
 
622
                                } else if (AtomSymbol.equals("F")) {
 
623
                                        factors[0] = 7.14/*7.34*/;
 
624
                                        factors[1] = 13.86;
 
625
                                        factors[2] = 5.68;
 
626
                                } else if (AtomSymbol.equals("Cl")) {
 
627
                                        factors[0] = 6.51;/*6.50*/
 
628
                                        factors[1] = 11.02;
 
629
                                        factors[2] = 4.52;
 
630
                                } else if (AtomSymbol.equals("Br")) {
 
631
                                        factors[0] = 5.20;
 
632
                                        factors[1] = 9.68;
 
633
                                        factors[2] = 4.48;
 
634
                                } else if (AtomSymbol.equals("I")) {
 
635
                                        factors[0] = 4.95;
 
636
                                        factors[1] = 8.81;
 
637
                                        factors[2] = 3.86;
 
638
                                }
 
639
                        
 
640
                                gasteigerFactors[k][STEP_SIZE * i + i] = factors[0];
 
641
                                gasteigerFactors[k][STEP_SIZE * i + i + 1] = factors[1];
 
642
                                gasteigerFactors[k][STEP_SIZE * i + i + 2] = factors[2];
 
643
                                gasteigerFactors[k][STEP_SIZE * i + i + 5] = ac.getAtom(i).getCharge();
 
644
                                
 
645
                                if (factors[0] == 0 && factors[1] == 0 && factors[2] == 0) {
 
646
                                        gasteigerFactors[k][STEP_SIZE * i + i + 3] = 1;
 
647
                                } else {
 
648
                                        gasteigerFactors[k][STEP_SIZE * i + i + 3] = factors[0] + factors[1] + factors[2];
 
649
                                }
 
650
                        }
 
651
                }
 
652
                
 
653
 
 
654
                return gasteigerFactors;
 
655
        }
 
656
        /**
 
657
     * clean the flags CDKConstants.REACTIVE_CENTER from the molecule
 
658
     * 
 
659
         * @param ac
 
660
         */
 
661
        public void cleanFlagReactiveCenter(IAtomContainer ac){
 
662
                for(int j = 0 ; j < ac.getAtomCount(); j++)
 
663
                        ac.getAtom(j).setFlag(CDKConstants.REACTIVE_CENTER, false);
 
664
                for(int j = 0 ; j < ac.getBondCount(); j++)
 
665
                        ac.getBond(j).setFlag(CDKConstants.REACTIVE_CENTER, false);
 
666
        }
 
667
}
 
668