ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/dipmom.c
Revision: 58
Committed: Mon Dec 1 06:44:59 2008 UTC (13 years, 8 months ago) by wdelano
File size: 4975 byte(s)
Log Message:
code merge 20081130
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "bonds_ff.h"
6
7 EXTERN struct t_dipolemom {
8 double total, xdipole, ydipole, zdipole;
9 } dipolemom;
10 struct t_quadmom {
11 double xx,xy,xz,yx,yy,yz,zx,zy,zz;
12 } quadmom;
13
14
15 void charge_dipole(void);
16 void dipole_dipole(void);
17 void vibcharge_dipole(double **);
18 void vibdipole_dipole(double **);
19
20 void charge_dipole(void)
21 {
22 int i;
23 double weight, xcenter, ycenter, zcenter;
24 double xsum, ysum, zsum, debye;
25
26 debye = 4.8033324;
27
28 weight = 0.0;
29 xcenter = 0.0;
30 ycenter = 0.0;
31 zcenter = 0.0;
32 for (i=1; i <= natom; i++)
33 {
34 weight += atom[i].atomwt;
35 xcenter += atom[i].x*atom[i].atomwt;
36 ycenter += atom[i].y*atom[i].atomwt;
37 zcenter += atom[i].z*atom[i].atomwt;
38 }
39 xcenter /= weight;
40 ycenter /= weight;
41 zcenter /= weight;
42
43 for(i=1; i <= natom; i++)
44 {
45 xsum += (atom[i].x-xcenter)*atom[i].charge;
46 ysum += (atom[i].y-ycenter)*atom[i].charge;
47 zsum += (atom[i].z-zcenter)*atom[i].charge;
48 }
49 dipolemom.xdipole = xsum*debye;
50 dipolemom.ydipole = ysum*debye;
51 dipolemom.zdipole = zsum*debye;
52 }
53
54 void dipole_dipole()
55 {
56 int i,ia, ib;
57 double xsum, ysum, zsum;
58 double xi, yi, zi, fik, ri2;
59
60 xsum = 0.0;
61 ysum = 0.0;
62 zsum = 0.0;
63 for (i=0; i < bonds_ff.nbnd ; i++)
64 {
65 if (bonds_ff.bmom[i] != 0.0)
66 {
67 ia = bonds_ff.idpl[i][0];
68 ib = bonds_ff.idpl[i][1];
69 xi = atom[ib].x - atom[ia].x;
70 yi = atom[ib].y - atom[ia].y;
71 zi = atom[ib].z - atom[ia].z;
72 ri2 = xi*xi + yi*yi + zi*zi;
73 fik = bonds_ff.bmom[i]/sqrt(ri2);
74 xsum += fik*xi;
75 ysum += fik*yi;
76 zsum += fik*zi;
77 }
78 }
79 dipolemom.xdipole = xsum;
80 dipolemom.ydipole = ysum;
81 dipolemom.zdipole = zsum;
82
83 }
84 /* ================================================= */
85 void vibcharge_dipole(double **acoord)
86 {
87 int i;
88 double weight, xcenter, ycenter, zcenter;
89 double xsum, ysum, zsum, debye;
90 double xy,xz,yx,yz,zx,zy;
91
92 debye = 4.8033324;
93
94 weight = 0.0;
95 xcenter = 0.0;
96 ycenter = 0.0;
97 zcenter = 0.0;
98 for (i=1; i <= natom; i++)
99 {
100 weight += atom[i].atomwt;
101 xcenter += acoord[i][0]*atom[i].atomwt;
102 ycenter += acoord[i][1]*atom[i].atomwt;
103 zcenter += acoord[i][2]*atom[i].atomwt;
104 }
105 xcenter /= weight;
106 ycenter /= weight;
107 zcenter /= weight;
108
109 xsum = 0.0;
110 ysum = 0.0;
111 zsum = 0.0;
112 for(i=1; i <= natom; i++)
113 {
114 xsum += (acoord[i][0]-xcenter)*atom[i].charge;
115 ysum += (acoord[i][1]-ycenter)*atom[i].charge;
116 zsum += (acoord[i][2]-zcenter)*atom[i].charge;
117 }
118 dipolemom.xdipole = xsum*debye;
119 dipolemom.ydipole = ysum*debye;
120 dipolemom.zdipole = zsum*debye;
121 // quadrapole moment
122 xsum = 0.0;
123 ysum = 0.0;
124 zsum = 0.0;
125 xy = xz = yx = yz = zx = zy = 0.0;
126 for(i=1; i <= natom; i++)
127 {
128 xsum += (acoord[i][0]-xcenter)*(acoord[i][0]-xcenter)*atom[i].charge;
129 ysum += (acoord[i][1]-ycenter)*(acoord[i][1]-ycenter)*atom[i].charge;
130 zsum += (acoord[i][2]-zcenter)*(acoord[i][2]-zcenter)*atom[i].charge;
131 xy += (acoord[i][0]-xcenter)*(acoord[i][1]-ycenter)*atom[i].charge;
132 xz += (acoord[i][0]-xcenter)*(acoord[i][2]-zcenter)*atom[i].charge;
133 yx += (acoord[i][1]-ycenter)*(acoord[i][0]-xcenter)*atom[i].charge;
134 yz += (acoord[i][1]-ycenter)*(acoord[i][2]-zcenter)*atom[i].charge;
135 zx += (acoord[i][2]-zcenter)*(acoord[i][0]-xcenter)*atom[i].charge;
136 zy += (acoord[i][2]-zcenter)*(acoord[i][1]-ycenter)*atom[i].charge;
137 }
138 quadmom.xx = xsum;
139 quadmom.xy = xy;
140 quadmom.xz = xz;
141 quadmom.yy = ysum;
142 quadmom.yx = yx;
143 quadmom.yz = yz;
144 quadmom.zz = zsum;
145 quadmom.zx = zx;
146 quadmom.zy = zy;
147
148 }
149
150 void vibdipole_dipole(double **acoord)
151 {
152 int i,ia, ib;
153 double xsum, ysum, zsum,bmomj;
154 double xi, yi, zi, fik, ri2;
155
156 xsum = 0.0;
157 ysum = 0.0;
158 zsum = 0.0;
159 for (i=0; i < bonds_ff.nbnd ; i++)
160 {
161 ia = bonds_ff.i12[i][0];
162 ib = bonds_ff.i12[i][1];
163 if (atom[ia].type == 1 && atom[ib].type == 5)
164 bmomj = -0.25;
165 else if (atom[ia].type == 1 && atom[ib].type == 36)
166 bmomj = -0.25;
167 else
168 bmomj = bonds_ff.bmom[i];
169 if (bmomj != 0.0)
170 {
171 xi = acoord[ib][0] - acoord[ia][0];
172 yi = acoord[ib][1] - acoord[ia][1];
173 zi = acoord[ib][2] - acoord[ia][2];
174 ri2 = xi*xi + yi*yi + zi*zi;
175 fik = bmomj/sqrt(ri2);
176 xsum += fik*xi;
177 ysum += fik*yi;
178 zsum += fik*zi;
179 }
180 }
181 dipolemom.xdipole = xsum;
182 dipolemom.ydipole = ysum;
183 dipolemom.zdipole = zsum;
184
185 }
186