← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/org/openscience/cdk/math/qm/ClosedShellJob.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/math/qm/ClosedShellJob.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) 2001-2007  The Chemistry Development Kit (CDK) project
 
7
 * 
 
8
 * Contact: cdk-devel@lists.sf.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
 */
 
29
package org.openscience.cdk.math.qm;
 
30
 
 
31
import org.openscience.cdk.math.Matrix;
 
32
import org.openscience.cdk.math.Vector;
 
33
import org.openscience.cdk.tools.LoggingTool;
 
34
 
 
35
/**
 
36
 * Calculates the orbitals and orbital energies of electron systems
 
37
 * with closed shells
 
38
 * 
 
39
 * @author  Stephan Michels <stephan@vern.chem.tu-berlin.de>
 
40
 * @cdk.created 2001-06-14
 
41
 */
 
42
public class ClosedShellJob
 
43
{
 
44
  private Orbitals orbitals;
 
45
  private Vector E;
 
46
 
 
47
  private LoggingTool log = new LoggingTool(ClosedShellJob.class);
 
48
 
 
49
  private int iterations = 0;
 
50
 
 
51
  public ClosedShellJob(Orbitals orbitals)
 
52
  { 
 
53
    this.orbitals = orbitals;
 
54
  }
 
55
 
 
56
  public Vector getEnergies()
 
57
  {
 
58
    return E.duplicate();
 
59
  }
 
60
 
 
61
  /**
 
62
   * Sorts the orbitals by their energies
 
63
   */
 
64
  private void sort(Matrix C, Vector E)
 
65
  {
 
66
    int i,j;
 
67
    double value;
 
68
    boolean changed;
 
69
    do
 
70
    {
 
71
      changed = false;
 
72
      for(i=1; i<E.size; i++)
 
73
        if (E.vector[i-1]>E.vector[i])
 
74
        {
 
75
          value = E.vector[i];
 
76
          E.vector[i] = E.vector[i-1];
 
77
          E.vector[i-1] = value;
 
78
 
 
79
          for(j=0; j<C.rows; j++)
 
80
          {
 
81
            value = C.matrix[j][i];
 
82
            C.matrix[j][i] = C.matrix[j][i-1];
 
83
            C.matrix[j][i-1] = value;
 
84
          }
 
85
          changed = true;
 
86
        }
 
87
    } while (changed);
 
88
  }
 
89
 
 
90
  private Matrix calculateS(IBasis basis)
 
91
  {
 
92
    int size = basis.getSize();
 
93
    Matrix S = new Matrix(size,size);
 
94
    int i,j;
 
95
    for(i=0; i<size; i++)
 
96
      for(j=0; j<size; j++)
 
97
        S.matrix[i][j] = basis.calcS(i,j);
 
98
 
 
99
    return S;
 
100
  }
 
101
 
 
102
  /**
 
103
   * Calculates the matrix for the kinetic energy
 
104
   *
 
105
   * T_i,j = (1/2) * -<d^2/dx^2 chi_i | chi_j>
 
106
   */
 
107
  private Matrix calculateT(IBasis basis)
 
108
  {
 
109
    int size = basis.getSize();
 
110
    Matrix J = new Matrix(size,size);
 
111
    int i,j; 
 
112
    for(i=0; i<size; i++)
 
113
      for(j=0; j<size; j++)
 
114
        // (1/2) * -<d^2/dx^2 chi_i | chi_j>
 
115
        J.matrix[i][j] = basis.calcJ(j,i)/2; // Vorsicht indizies sind vertauscht
 
116
 
 
117
    return J;
 
118
  }
 
119
 
 
120
  /**
 
121
   * Calculates the matrix for the potential matrix
 
122
   *
 
123
   * V_i,j = <chi_i | 1/r | chi_j>
 
124
   */
 
125
  private Matrix calculateV(IBasis basis)
 
126
  {
 
127
    int size = basis.getSize();
 
128
    Matrix V = new Matrix(size,size);
 
129
    int i,j;
 
130
    for(i=0; i<size; i++)
 
131
      for(j=0; j<size; j++)
 
132
        V.matrix[i][j] = basis.calcV(i,j);
 
133
 
 
134
    return V;
 
135
  }
 
136
 
 
137
  /**
 
138
   * Calculates thes values for the 2 electron interactions
 
139
   */
 
140
  private double[][][][] calculateI(IBasis basis)
 
141
  {
 
142
    int i,j,k,l;
 
143
    int size = basis.getSize();
 
144
 
 
145
    double[][][][] result = new double[size][][][];
 
146
    for(i=0; i<size; i++)
 
147
    {
 
148
      result[i] = new double[i+1][][];
 
149
      for(j=0; j<=i; j++)
 
150
      {
 
151
        result[i][j] = new double[size][];
 
152
        for(k=0; k<size; k++) 
 
153
        {
 
154
          result[i][j][k] = new double[k+1];
 
155
          for(l=0; l<=k; l++)
 
156
          {
 
157
            result[i][j][k][l] = basis.calcI(i,j,k,l);
 
158
            //log.println("("+(i+1)+" "+(j+1)+"|"+(k+1)+" "+(l+1)+")="+result[i][j][k][l]);
 
159
          }
 
160
        }
 
161
      }
 
162
    }
 
163
    return result;
 
164
  }
 
165
 
 
166
  /**
 
167
   * Calculates the density matrix
 
168
   */
 
169
  private Matrix calculateD(IBasis basis, Matrix C, int count_electrons)
 
170
  {
 
171
    int i,j,k;
 
172
    int size = basis.getSize();
 
173
    int orbitals = C.getColumns();
 
174
    int occ = count_electrons/2;
 
175
    int locc = count_electrons%2;
 
176
    Matrix D = new Matrix(size,size);
 
177
    log.debug("D:occ="+occ+" locc="+locc);
 
178
 
 
179
//    if (locc!=0)
 
180
//      logger.debug("This class work only correct for closed shells");
 
181
 
 
182
    for(i=0; i<size; i++)
 
183
      for(j=0; j<size; j++)
 
184
      {
 
185
        D.matrix[i][j] = 0d;
 
186
        for(k=0; (k<orbitals) && (k<occ); k++)
 
187
          D.matrix[i][j] += 2d*C.matrix[i][k]*C.matrix[j][k];
 
188
 
 
189
        if ((locc==1) && (k+1<orbitals))
 
190
          D.matrix[i][j] += C.matrix[i][k+1]*C.matrix[j][k+1];
 
191
      }
 
192
    return D;
 
193
  }
 
194
 
 
195
  private Matrix calculateJ(IBasis basis, double[][][][] I, Matrix D)
 
196
  { 
 
197
    int i,j,k,l;
 
198
    int size = basis.getSize();
 
199
    Matrix J = new Matrix(size,size);
 
200
    for(i=0; i<size; i++)
 
201
      for(j=0; j<size; j++)
 
202
      {
 
203
        J.matrix[i][j] = 0;
 
204
        for(k=0; k<size; k++)
 
205
          for(l=0; l<size; l++)
 
206
          {
 
207
            if (i>=j)
 
208
            {
 
209
              if (k>=l)
 
210
                J.matrix[i][j] += D.matrix[k][l]*I[i][j][k][l];
 
211
              else
 
212
                J.matrix[i][j] += D.matrix[k][l]*I[i][j][l][k];
 
213
            }
 
214
            else
 
215
            {
 
216
              if (k>=l)
 
217
                J.matrix[i][j] += D.matrix[k][l]*I[j][i][k][l];
 
218
              else
 
219
                J.matrix[i][j] += D.matrix[k][l]*I[j][i][l][k];
 
220
            }
 
221
          }
 
222
        J.matrix[i][j] *= 2d;
 
223
      }
 
224
    return J;
 
225
  }
 
226
 
 
227
  private Matrix calculateK(IBasis basis, double[][][][] I, Matrix D)
 
228
  { 
 
229
    int i,j,k,l;
 
230
    int size = basis.getSize();
 
231
    Matrix K = new Matrix(size,size);
 
232
    for(i=0; i<size; i++)
 
233
      for(j=0; j<size; j++)
 
234
      { 
 
235
        K.matrix[i][j] = 0;
 
236
        for(k=0; k<size; k++)
 
237
          for(l=0; l<size; l++)
 
238
          { 
 
239
            if (i>=j)
 
240
            {
 
241
              if (k>=l)
 
242
                K.matrix[i][j] += D.matrix[k][l]*I[i][j][k][l];
 
243
              else
 
244
                K.matrix[i][j] += D.matrix[k][l]*I[i][j][l][k];
 
245
            }
 
246
            else
 
247
            {
 
248
              if (k>=l)
 
249
                K.matrix[i][j] += D.matrix[k][l]*I[j][i][k][l];
 
250
              else
 
251
                K.matrix[i][j] += D.matrix[k][l]*I[j][i][l][k];
 
252
            }
 
253
          } 
 
254
      }   
 
255
    return K;
 
256
  }
 
257
 
 
258
  private double contraction(Matrix A, Matrix B)
 
259
  {
 
260
    int i,j;
 
261
    double result = 0;
 
262
    for(i=0; i<A.rows; i++)
 
263
      for(j=0; j<A.columns; j++)
 
264
        result += A.matrix[i][j]*B.matrix[i][j];
 
265
    return result;
 
266
  }
 
267
  
 
268
  public Orbitals calculate()
 
269
  {
 
270
    long time = System.currentTimeMillis();
 
271
 
 
272
    Matrix C,S,T,V,HAO,H,D,J,K,F,U;
 
273
    double[][][][] I;
 
274
    double energy;
 
275
    IBasis basis = orbitals.getBasis();
 
276
 
 
277
    int count_electrons = orbitals.getCountElectrons();
 
278
 
 
279
    C = orbitals.getCoefficients().duplicate();
 
280
 
 
281
    S = calculateS(basis);
 
282
 
 
283
    log.debug("S = \n"+S+"\n");
 
284
 
 
285
    log.debug("C = \n"+C+"\n");
 
286
 
 
287
    C = C.orthonormalize(S);
 
288
    log.debug("C' = \n"+C+"\n");
 
289
    log.debug("C't * S * C' = \n"+S.similar(C)+"\n");
 
290
 
 
291
    T = calculateT(basis);
 
292
    log.debug("T = \n"+T+"\n");
 
293
 
 
294
    V = calculateV(basis);
 
295
    log.debug("V = \n"+V+"\n");
 
296
 
 
297
    HAO = T.add(V);
 
298
    log.debug("HAO = \n"+HAO+"\n");
 
299
 
 
300
    H = HAO.similar(C);
 
301
    log.debug("H = C't * HAO * C' = \n"+H.similar(C)+"\n");
 
302
 
 
303
    U = H.diagonalize(50);
 
304
    E = H.similar(U).getVectorFromDiagonal();
 
305
    C = C.mul(U);
 
306
    sort(C,E);
 
307
    log.debug("C(neu) = \n"+C+"\n");
 
308
    log.debug("E = \n"+E+"\n");
 
309
 
 
310
    for(int j=0; j<E.size; j++)
 
311
      log.debug("E("+(j+1)+".Orbital)="+(E.vector[j]*27.211)+" eV");
 
312
 
 
313
    time = System.currentTimeMillis()-time;
 
314
    log.debug("Time = "+time+" ms");
 
315
    time = System.currentTimeMillis();
 
316
 
 
317
    if (iterations>0)
 
318
      I = calculateI(basis);
 
319
    else
 
320
      I = null;
 
321
 
 
322
    for(int i=0; i<iterations; i++)
 
323
    {
 
324
      log.debug((i+1)+".Durchlauf\n");
 
325
 
 
326
      time = System.currentTimeMillis();
 
327
 
 
328
      log.debug("C't * S * C' = \n"+S.similar(C)+"\n");
 
329
 
 
330
      log.debug("count of electrons = "+count_electrons+"\n");
 
331
 
 
332
      D = calculateD(basis, C, count_electrons);
 
333
      log.debug("D = \n"+D+"\n");
 
334
 
 
335
      //log.println("2*contraction(D*S) = "+(D.mul(S)).contraction()*2+"\n");
 
336
      log.debug("2*contraction(D*S) = "+contraction(D,S)*2+"\n");
 
337
 
 
338
      //J = calculateJ(basis, D);
 
339
      J = calculateJ(basis, I, D);
 
340
      log.debug("J = \n"+J+"\n");
 
341
 
 
342
      //K = calculateK(basis, D);
 
343
      K = calculateK(basis, I, D);
 
344
      log.debug("K = \n"+K+"\n");
 
345
 
 
346
      F = HAO.add(J).sub(K);
 
347
      log.debug("F = H+J-K = \n"+F+"\n");
 
348
      
 
349
      H = F.similar(C);
 
350
      log.debug("H = C't * F * C' = \n"+H+"\n");
 
351
 
 
352
      U = H.diagonalize(50);
 
353
      E = H.similar(U).getVectorFromDiagonal();
 
354
      C = C.mul(U);
 
355
      sort(C,E);
 
356
      log.debug("C(neu) = \n"+C+"\n");
 
357
      log.debug("E = \n"+E+"\n");
 
358
 
 
359
      for(int j=0; j<E.size; j++)
 
360
        log.debug("E("+(j+1)+".Orbital)="+(E.vector[j]*27.211)+" eV");
 
361
 
 
362
      energy = contraction(D,HAO.add(F));
 
363
      log.debug("Gesamtenergie = "+energy+" ("+energy*27.211+" eV)\n");
 
364
 
 
365
      time = System.currentTimeMillis()-time;
 
366
      log.debug("Time = "+time+" ms");
 
367
 
 
368
      System.gc();
 
369
    }
 
370
 
 
371
    return new Orbitals(basis, C);
 
372
  }
 
373
}