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 (12 years, 4 months ago) by tjod
File size: 5587 byte(s)
Log Message:
test

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    
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