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