ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/dipmom.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 5587 byte(s)
Log Message:
test

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 EXTERN struct t_dmomv {
11 float xn,yn,zn,xp,yp,zp;
12 } dmomv;
13 struct t_quadmom {
14 double xx,xy,xz,yx,yy,yz,zx,zy,zz;
15 } quadmom;
16
17
18 void charge_dipole(void);
19 void dipole_dipole(void);
20 void vibcharge_dipole(double **);
21 void vibdipole_dipole(double **);
22
23 void charge_dipole(void)
24 {
25 int i;
26 double weight, xcenter, ycenter, zcenter;
27 double xsum, ysum, zsum, debye;
28 double cu;
29
30 debye = 4.8033324;
31
32 weight = 0.0;
33 xcenter = 0.0;
34 ycenter = 0.0;
35 zcenter = 0.0;
36 for (i=1; i <= natom; i++)
37 {
38 weight += atom[i].atomwt;
39 xcenter += atom[i].x*atom[i].atomwt;
40 ycenter += atom[i].y*atom[i].atomwt;
41 zcenter += atom[i].z*atom[i].atomwt;
42 }
43 xcenter /= weight;
44 ycenter /= weight;
45 zcenter /= weight;
46
47 xsum = 0.0;
48 ysum = 0.0;
49 zsum = 0.0;
50 dmomv.xp = 0.0;
51 dmomv.yp = 0.0;
52 dmomv.xp = 0.0;
53 dmomv.xn = 0.0;
54 dmomv.yn = 0.0;
55 dmomv.zn = 0.0;
56 for(i=1; i <= natom; i++)
57 {
58 xsum += (atom[i].x-xcenter)*atom[i].charge;
59 ysum += (atom[i].y-ycenter)*atom[i].charge;
60 zsum += (atom[i].z-zcenter)*atom[i].charge;
61
62 cu = atom[i].charge;
63 if ( cu < 0.0)
64 {
65 dmomv.xn -= atom[i].x*cu;
66 dmomv.yn -= atom[i].y*cu;
67 dmomv.zn -= atom[i].z*cu;
68 } else if (cu > 0.0)
69 {
70 dmomv.xp += atom[i].x*cu;
71 dmomv.yp += atom[i].y*cu;
72 dmomv.zp += atom[i].z*cu;
73 }
74
75 }
76 dipolemom.xdipole = xsum*debye;
77 dipolemom.ydipole = ysum*debye;
78 dipolemom.zdipole = zsum*debye;
79 }
80
81 void dipole_dipole()
82 {
83 int i,ia, ib;
84 double xsum, ysum, zsum;
85 double xi, yi, zi, fik, ri2;
86
87 xsum = 0.0;
88 ysum = 0.0;
89 zsum = 0.0;
90 for (i=0; i < bonds_ff.nbnd ; i++)
91 {
92 if (bonds_ff.bmom[i] != 0.0)
93 {
94 ia = bonds_ff.idpl[i][0];
95 ib = bonds_ff.idpl[i][1];
96 xi = atom[ib].x - atom[ia].x;
97 yi = atom[ib].y - atom[ia].y;
98 zi = atom[ib].z - atom[ia].z;
99 ri2 = xi*xi + yi*yi + zi*zi;
100 fik = bonds_ff.bmom[i]/sqrt(ri2);
101 xsum += fik*xi;
102 ysum += fik*yi;
103 zsum += fik*zi;
104 }
105 }
106 dipolemom.xdipole = xsum;
107 dipolemom.ydipole = ysum;
108 dipolemom.zdipole = zsum;
109
110 }
111 /* ================================================= */
112 void vibcharge_dipole(double **acoord)
113 {
114 int i;
115 double weight, xcenter, ycenter, zcenter;
116 double xsum, ysum, zsum, debye;
117 double xy,xz,yx,yz,zx,zy;
118
119 debye = 4.8033324;
120
121 weight = 0.0;
122 xcenter = 0.0;
123 ycenter = 0.0;
124 zcenter = 0.0;
125 for (i=1; i <= natom; i++)
126 {
127 weight += atom[i].atomwt;
128 xcenter += acoord[i][0]*atom[i].atomwt;
129 ycenter += acoord[i][1]*atom[i].atomwt;
130 zcenter += acoord[i][2]*atom[i].atomwt;
131 }
132 xcenter /= weight;
133 ycenter /= weight;
134 zcenter /= weight;
135
136 xsum = 0.0;
137 ysum = 0.0;
138 zsum = 0.0;
139 for(i=1; i <= natom; i++)
140 {
141 xsum += (acoord[i][0]-xcenter)*atom[i].charge;
142 ysum += (acoord[i][1]-ycenter)*atom[i].charge;
143 zsum += (acoord[i][2]-zcenter)*atom[i].charge;
144 }
145 dipolemom.xdipole = xsum*debye;
146 dipolemom.ydipole = ysum*debye;
147 dipolemom.zdipole = zsum*debye;
148 // quadrapole moment
149 xsum = 0.0;
150 ysum = 0.0;
151 zsum = 0.0;
152 xy = xz = yx = yz = zx = zy = 0.0;
153 for(i=1; i <= natom; i++)
154 {
155 xsum += (acoord[i][0]-xcenter)*(acoord[i][0]-xcenter)*atom[i].charge;
156 ysum += (acoord[i][1]-ycenter)*(acoord[i][1]-ycenter)*atom[i].charge;
157 zsum += (acoord[i][2]-zcenter)*(acoord[i][2]-zcenter)*atom[i].charge;
158 xy += (acoord[i][0]-xcenter)*(acoord[i][1]-ycenter)*atom[i].charge;
159 xz += (acoord[i][0]-xcenter)*(acoord[i][2]-zcenter)*atom[i].charge;
160 yx += (acoord[i][1]-ycenter)*(acoord[i][0]-xcenter)*atom[i].charge;
161 yz += (acoord[i][1]-ycenter)*(acoord[i][2]-zcenter)*atom[i].charge;
162 zx += (acoord[i][2]-zcenter)*(acoord[i][0]-xcenter)*atom[i].charge;
163 zy += (acoord[i][2]-zcenter)*(acoord[i][1]-ycenter)*atom[i].charge;
164 }
165 quadmom.xx = xsum;
166 quadmom.xy = xy;
167 quadmom.xz = xz;
168 quadmom.yy = ysum;
169 quadmom.yx = yx;
170 quadmom.yz = yz;
171 quadmom.zz = zsum;
172 quadmom.zx = zx;
173 quadmom.zy = zy;
174
175 }
176
177 void vibdipole_dipole(double **acoord)
178 {
179 int i,ia, ib;
180 double xsum, ysum, zsum,bmomj;
181 double xi, yi, zi, fik, ri2;
182
183 xsum = 0.0;
184 ysum = 0.0;
185 zsum = 0.0;
186 for (i=0; i < bonds_ff.nbnd ; i++)
187 {
188 ia = bonds_ff.i12[i][0];
189 ib = bonds_ff.i12[i][1];
190 if (atom[ia].type == 1 && atom[ib].type == 5)
191 bmomj = -0.25;
192 else if (atom[ia].type == 1 && atom[ib].type == 36)
193 bmomj = -0.25;
194 else
195 bmomj = bonds_ff.bmom[i];
196 if (bmomj != 0.0)
197 {
198 xi = acoord[ib][0] - acoord[ia][0];
199 yi = acoord[ib][1] - acoord[ia][1];
200 zi = acoord[ib][2] - acoord[ia][2];
201 ri2 = xi*xi + yi*yi + zi*zi;
202 fik = bmomj/sqrt(ri2);
203 xsum += fik*xi;
204 ysum += fik*yi;
205 zsum += fik*zi;
206 }
207 }
208 dipolemom.xdipole = xsum;
209 dipolemom.ydipole = ysum;
210 dipolemom.zdipole = zsum;
211
212 }
213