← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cints/DFT/weighting.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cints/DFT/weighting.cc
 
1
/*! \file weighting.cc
 
2
    \ingroup CINTS
 
3
    \author Shawn Brown
 
4
   
 
5
    This code contains some functions used for 
 
6
    the calculation of the weighting functions
 
7
    for the Becke scheme.
 
8
    
 
9
    Ref.  Becke, J. Chem. Phys., Vol. 88, pg. 2547, 1988.
 
10
    
 
11
    ----------------------------------------------*/
 
12
#include <cmath>
 
13
#include <cstdio>
 
14
#include <memory.h>
 
15
#include <cstdlib>
 
16
#include <libipv1/ip_lib.h>
 
17
#include <libciomr/libciomr.h>
 
18
#include <libint/libint.h>
 
19
 
 
20
#include"defines.h"
 
21
#define EXTERN
 
22
#include"global.h"
 
23
#include <stdexcept>
 
24
#include"small_fns.h"
 
25
 
 
26
/* Declare functions in this code */
 
27
 
 
28
namespace psi { namespace CINTS {
 
29
double u_calc(int i, int j, struct coordinates geom); 
 
30
double v_calc(int i, int j, double uij);
 
31
double f_u(double vij);
 
32
double s_u(double f);
 
33
 
 
34
double weight_calc(int atomn,struct coordinates  geom,int k_order){
 
35
    int i,j,k,l;
 
36
    int natoms;
 
37
    double utemp,vtemp,ftemp;
 
38
    double **s_mat;
 
39
    double *ptemp;
 
40
    double sum=0.0;
 
41
    double weight;
 
42
    
 
43
    natoms = Molecule.num_atoms;
 
44
    /*natoms = Symmetry.num_unique_atoms;*/
 
45
    
 
46
    ptemp = init_array(natoms);
 
47
    s_mat = block_matrix(natoms,natoms);
 
48
    
 
49
    /* calculate all possible s's*/
 
50
    for(i=0;i<natoms;i++){
 
51
        for(j=0;j<natoms;j++){
 
52
            if(i!=j){
 
53
                utemp = u_calc(i,j,geom);
 
54
                vtemp = v_calc(i,j,utemp);
 
55
                ftemp = f_u(utemp);
 
56
                for(k=1;k<k_order;k++){
 
57
                    ftemp = f_u(ftemp);
 
58
                }
 
59
                s_mat[i][j] = s_u(ftemp);
 
60
            }
 
61
        }
 
62
    }
 
63
        
 
64
    /* calculate all of the cell functions into an array */
 
65
    
 
66
    for(i=0;i<natoms;i++){
 
67
        ptemp[i] = 1.0;
 
68
        for(j=0;j<natoms;j++){
 
69
            if(i!=j)
 
70
                ptemp[i] *= s_mat[i][j];
 
71
        }
 
72
        sum += ptemp[i];
 
73
    }
 
74
    
 
75
    weight=ptemp[atomn]/sum;
 
76
    free(ptemp);
 
77
    free_block(s_mat);
 
78
    
 
79
    return weight;
 
80
 
 
81
}
 
82
                
 
83
double u_calc(int i, int j, struct coordinates geom){
 
84
    
 
85
    double ri;
 
86
    double rj;
 
87
    double rij;
 
88
    /* calculate the distance of the grid point from
 
89
       each atom */
 
90
    
 
91
    ri = distance_calc(Molecule.centers[i],geom);
 
92
    rj = distance_calc(Molecule.centers[j],geom);
 
93
    rij = distance_calc(Molecule.centers[i],Molecule.centers[j]);
 
94
    return (ri-rj)/rij;
 
95
}
 
96
 
 
97
 
 
98
double f_u(double vij){
 
99
    
 
100
    return 1.5*vij-0.5*(vij*vij*vij);
 
101
}
 
102
 
 
103
double s_u(double f){
 
104
    
 
105
    return 0.5*(1.0-f);
 
106
}
 
107
    
 
108
double v_calc(int atomi, int atomj,double uij){
 
109
    
 
110
    double aij;
 
111
    double utmp;
 
112
    double rattmp;
 
113
    double tmp1,tmp2,tmp3;
 
114
 
 
115
    rattmp = DFT_options.grid.atomic_grid[atomi].Bragg_radii/
 
116
        DFT_options.grid.atomic_grid[atomj].Bragg_radii;
 
117
    utmp = (rattmp-1)/(rattmp+1);
 
118
    aij=utmp/((utmp*utmp)-1);
 
119
    if(aij>0.5) aij = 0.5;
 
120
    if(aij<-0.5) aij = -0.5;
 
121
    tmp1 = uij*uij;
 
122
    tmp2 = aij*(1-tmp1);
 
123
    tmp3 = uij+tmp2;
 
124
 
 
125
    return uij+tmp2;
 
126
}               
 
127
};}