← Back to branch summary

~ubuntu-branches/ubuntu/trusty/cdk/trusty-proposed

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/geometry/RDFCalculator.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/geometry/RDFCalculator.java
 
1
/* $RCSfile$
 
2
 * $Author: egonw $    
 
3
 * $Date: 2007-01-04 18:46:10 +0100 (Thu, 04 Jan 2007) $    
 
4
 * $Revision: 7636 $
 
5
 * 
 
6
 * Copyright (C) 2005-2007  The Chemistry Development Kit (CDK) project
 
7
 * 
 
8
 * Contact: cdk-devel@lists.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
 * All we ask is that proper credit is given for our work, which includes
 
15
 * - but is not limited to - adding the above copyright notice to the beginning
 
16
 * of your source code files, and to any copyright notice that you may distribute
 
17
 * with programs based on this work.
 
18
 * 
 
19
 * This program is distributed in the hope that it will be useful,
 
20
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
21
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
22
 * GNU Lesser General Public License for more details.
 
23
 * 
 
24
 * You should have received a copy of the GNU Lesser General Public License
 
25
 * along with this program; if not, write to the Free Software
 
26
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
 
27
 */
 
28
package org.openscience.cdk.geometry;
 
29
 
 
30
import javax.vecmath.Point3d;
 
31
 
 
32
import org.openscience.cdk.Atom;
 
33
import org.openscience.cdk.AtomContainer;
 
34
import org.openscience.cdk.tools.LoggingTool;
 
35
import org.openscience.cdk.interfaces.IAtom;
 
36
 
 
37
/**
 
38
 * Calculator of radial distribution functions. The RDF has bins defined around
 
39
 * a point, i.e. the first bin starts at 0 Å and ends at 0.5*resolution
 
40
 * Å, and the second bins ends at 1.5*resulution Å.
 
41
 *
 
42
 * <p>By default, the RDF is unweighted. By implementing and registering a
 
43
 * <code>RDFWeightFunction</code>, the RDF can become weighted. For example,
 
44
 * to weight according to partial charge interaction, this code could be used:
 
45
 * <pre>
 
46
 * RDFCalculator calculator = new RDFCalculator(0.0, 5.0, 0.1, 0.0,
 
47
 *     new RDFWeightFunction() {
 
48
 *         public double calculate(Atom atom, Atom atom2) {
 
49
 *             return atom.getCharge()*atom2.getCharge();
 
50
 *         }
 
51
 *     }
 
52
 * );
 
53
 * </pre>
 
54
 *
 
55
 * @cdk.module  extra
 
56
 *
 
57
 * @author      Egon Willighagen
 
58
 * @cdk.created 2005-01-10
 
59
 *
 
60
 * @cdk.keyword radial distribution function
 
61
 * @cdk.keyword RDF
 
62
 *
 
63
 * @see         org.openscience.cdk.geometry.IRDFWeightFunction
 
64
 */
 
65
public class RDFCalculator {
 
66
 
 
67
    private LoggingTool logger;
 
68
    
 
69
    private double startCutoff;
 
70
    private double cutoff;
 
71
    private double resolution;
 
72
    private double peakWidth;
 
73
    
 
74
    private IRDFWeightFunction weightFunction;
 
75
    
 
76
    /**
 
77
     * Constructs a RDF calculator that calculates a unweighted, digitized
 
78
     * RDF function.
 
79
     *
 
80
     * @param startCutoff radial length in &Aring;ngstrom at which the RDF starts
 
81
     * @param cutoff      radial length in &Aring;ngstrom at which the RDF stops
 
82
     * @param resolution  width of the bins
 
83
     * @param peakWidth   width of the gaussian applied to the peaks in &Aring;ngstrom
 
84
     */
 
85
    public RDFCalculator(double startCutoff, double cutoff, double resolution, 
 
86
                         double peakWidth) {
 
87
        this(startCutoff, cutoff, resolution, peakWidth, null);
 
88
    }
 
89
 
 
90
    /**
 
91
     * Constructs a RDF calculator that calculates a digitized
 
92
     * RDF function.
 
93
     *
 
94
     * @param startCutoff    radial length in &Aring;ngstrom at which the RDF starts
 
95
     * @param cutoff         radial length in &Aring;ngstrom at which the RDF stops
 
96
     * @param resolution     width of the bins
 
97
     * @param peakWidth      width of the gaussian applied to the peaks in &Aring;ngstrom
 
98
     * @param weightFunction the weight function. If null, then an unweighted RDF is
 
99
     *                       calculated
 
100
     */
 
101
    public RDFCalculator(double startCutoff, double cutoff, double resolution, 
 
102
                         double peakWidth, IRDFWeightFunction weightFunction) {
 
103
        logger = new LoggingTool(this);
 
104
        
 
105
         this.startCutoff = startCutoff;
 
106
         this.cutoff = cutoff;
 
107
         this.resolution = resolution;
 
108
         this.peakWidth = peakWidth;
 
109
         this.weightFunction = weightFunction;
 
110
    }
 
111
    
 
112
    /**
 
113
     * Calculates a RDF for <code>Atom</code> atom in the environment
 
114
     * of the atoms in the <code>AtomContainer</code>.
 
115
     */
 
116
    public double[] calculate(AtomContainer container, Atom atom) {
 
117
        int length = (int)((cutoff-startCutoff)/resolution) + 1;
 
118
        logger.debug("Creating RDF of length ", length);
 
119
 
 
120
        // the next we need for Gaussian smoothing
 
121
        int binsToFillOnEachSide = (int)(peakWidth*3.0/resolution);
 
122
        double sigmaSquare = Math.pow(peakWidth, 2.0);
 
123
        double[] factors = new double[binsToFillOnEachSide];
 
124
        for (int binCounter=0; binCounter<binsToFillOnEachSide; binCounter++) {
 
125
            factors[binCounter] = Math.exp(-1.0*(Math.pow(((double)binCounter)*resolution, 2.0))/sigmaSquare);
 
126
        }
 
127
        
 
128
        // this we need always
 
129
        double[] rdf = new double[length];
 
130
        double distance = 0.0;
 
131
        int index = 0;
 
132
        
 
133
        Point3d atomPoint = atom.getPoint3d();
 
134
        java.util.Iterator atomsInContainer = container.atoms();
 
135
        while (atomsInContainer.hasNext()) {
 
136
                IAtom atomInContainer = (IAtom)atomsInContainer.next();
 
137
            distance = atomPoint.distance(atomInContainer.getPoint3d());
 
138
            index = (int)((distance-startCutoff)/this.resolution);
 
139
            double weight = 1.0;
 
140
            if (weightFunction != null) {
 
141
                weight = weightFunction.calculate(atom, atomInContainer);
 
142
            }
 
143
            rdf[index] += weight; // unweighted
 
144
            if (this.peakWidth > 0.0) {
 
145
                // apply Gaussian smoothing
 
146
                for (int binCounter=1; binCounter<=binsToFillOnEachSide; binCounter++) {
 
147
                    if ((index - binCounter) >= 0) {
 
148
                        rdf[index - binCounter] += weight*factors[binCounter];
 
149
                    }
 
150
                    if ((index + binCounter) < length) {
 
151
                        rdf[index + binCounter] += weight*factors[binCounter];
 
152
                    }
 
153
                }
 
154
            }
 
155
        }
 
156
        return rdf;
 
157
    }
 
158
    
 
159
}
 
160
 
 
161
 
 
162