ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/dipmom.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (13 years, 8 months ago) by gilbertke
File size: 3695 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "bonds_ff.h"
6 gilbertke 63 #include "dipmom.h"
7 tjod 3
8    
9 gilbertke 63 static void vibcharge_dipole(double **);
10     static void vibdipole_dipole(double **);
11 tjod 3
12     void charge_dipole(void)
13     {
14     int i;
15     double weight, xcenter, ycenter, zcenter;
16     double xsum, ysum, zsum, debye;
17    
18     debye = 4.8033324;
19    
20     weight = 0.0;
21     xcenter = 0.0;
22     ycenter = 0.0;
23     zcenter = 0.0;
24     for (i=1; i <= natom; i++)
25     {
26     weight += atom[i].atomwt;
27     xcenter += atom[i].x*atom[i].atomwt;
28     ycenter += atom[i].y*atom[i].atomwt;
29     zcenter += atom[i].z*atom[i].atomwt;
30     }
31     xcenter /= weight;
32     ycenter /= weight;
33     zcenter /= weight;
34    
35     for(i=1; i <= natom; i++)
36     {
37     xsum += (atom[i].x-xcenter)*atom[i].charge;
38     ysum += (atom[i].y-ycenter)*atom[i].charge;
39     zsum += (atom[i].z-zcenter)*atom[i].charge;
40     }
41     dipolemom.xdipole = xsum*debye;
42     dipolemom.ydipole = ysum*debye;
43     dipolemom.zdipole = zsum*debye;
44     }
45    
46     void dipole_dipole()
47     {
48     int i,ia, ib;
49     double xsum, ysum, zsum;
50     double xi, yi, zi, fik, ri2;
51    
52     xsum = 0.0;
53     ysum = 0.0;
54     zsum = 0.0;
55     for (i=0; i < bonds_ff.nbnd ; i++)
56     {
57     if (bonds_ff.bmom[i] != 0.0)
58     {
59     ia = bonds_ff.idpl[i][0];
60     ib = bonds_ff.idpl[i][1];
61     xi = atom[ib].x - atom[ia].x;
62     yi = atom[ib].y - atom[ia].y;
63     zi = atom[ib].z - atom[ia].z;
64     ri2 = xi*xi + yi*yi + zi*zi;
65     fik = bonds_ff.bmom[i]/sqrt(ri2);
66     xsum += fik*xi;
67     ysum += fik*yi;
68     zsum += fik*zi;
69     }
70     }
71     dipolemom.xdipole = xsum;
72     dipolemom.ydipole = ysum;
73     dipolemom.zdipole = zsum;
74    
75     }
76     /* ================================================= */
77     void vibcharge_dipole(double **acoord)
78     {
79     int i;
80     double weight, xcenter, ycenter, zcenter;
81     double xsum, ysum, zsum, debye;
82     double xy,xz,yx,yz,zx,zy;
83    
84     debye = 4.8033324;
85    
86     weight = 0.0;
87     xcenter = 0.0;
88     ycenter = 0.0;
89     zcenter = 0.0;
90     for (i=1; i <= natom; i++)
91     {
92     weight += atom[i].atomwt;
93     xcenter += acoord[i][0]*atom[i].atomwt;
94     ycenter += acoord[i][1]*atom[i].atomwt;
95     zcenter += acoord[i][2]*atom[i].atomwt;
96     }
97     xcenter /= weight;
98     ycenter /= weight;
99     zcenter /= weight;
100    
101     xsum = 0.0;
102     ysum = 0.0;
103     zsum = 0.0;
104     for(i=1; i <= natom; i++)
105     {
106     xsum += (acoord[i][0]-xcenter)*atom[i].charge;
107     ysum += (acoord[i][1]-ycenter)*atom[i].charge;
108     zsum += (acoord[i][2]-zcenter)*atom[i].charge;
109     }
110     dipolemom.xdipole = xsum*debye;
111     dipolemom.ydipole = ysum*debye;
112 gilbertke 63 dipolemom.zdipole = zsum*debye;
113 tjod 3 }
114    
115     void vibdipole_dipole(double **acoord)
116     {
117     int i,ia, ib;
118     double xsum, ysum, zsum,bmomj;
119     double xi, yi, zi, fik, ri2;
120    
121     xsum = 0.0;
122     ysum = 0.0;
123     zsum = 0.0;
124     for (i=0; i < bonds_ff.nbnd ; i++)
125     {
126     ia = bonds_ff.i12[i][0];
127     ib = bonds_ff.i12[i][1];
128     if (atom[ia].type == 1 && atom[ib].type == 5)
129     bmomj = -0.25;
130     else if (atom[ia].type == 1 && atom[ib].type == 36)
131     bmomj = -0.25;
132     else
133     bmomj = bonds_ff.bmom[i];
134     if (bmomj != 0.0)
135     {
136     xi = acoord[ib][0] - acoord[ia][0];
137     yi = acoord[ib][1] - acoord[ia][1];
138     zi = acoord[ib][2] - acoord[ia][2];
139     ri2 = xi*xi + yi*yi + zi*zi;
140     fik = bmomj/sqrt(ri2);
141     xsum += fik*xi;
142     ysum += fik*yi;
143     zsum += fik*zi;
144     }
145     }
146     dipolemom.xdipole = xsum;
147     dipolemom.ydipole = ysum;
148     dipolemom.zdipole = zsum;
149    
150     }
151