ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/hessian.c
Revision: 58
Committed: Mon Dec 1 06:44:59 2008 UTC (12 years, 6 months ago) by wdelano
File size: 4808 byte(s)
Log Message:
code merge 20081130
Line File contents
1 #define EXTERN extern
2
3
4 #include "pcwin.h"
5 #include "pcmod.h"
6
7 #include "energies.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "utility.h"
11 #include "pot.h"
12
13
14 EXTERN double hesscut;
15
16 void piseq(int, int);
17 void hessian(int, double *,int *, int *, int *,double *);
18 void ebond2(int);
19 void eangle2(int);
20 void estrbnd2(int);
21 void eangang2(int);
22 void eopbend2(int);
23 void etorsion2(int);
24 void estrtor2(int);
25 void eimprop2(int);
26 void eimptors2(int);
27 void ehbond2(int);
28 void ebuck2(int);
29 void ecoord_bond2(int);
30 void elj2(int);
31 void echarge2(int);
32 void edipole2(int);
33 void ehal2(int);
34 void ebufcharge2(int);
35 void eopbend_wilson2(int);
36 void ewald2(int);
37 void eurey2(int);
38 void elj_qq2(int);
39 void esolv2(int);
40 void ehigh_coord2(int);
41 void egeom2(int);
42 void message_alert(char *,char *);
43
44 void hessian(int maxhess, double *h, int *hinit, int *hstop,int *hindex, double *hdiag)
45 {
46 int i,j,k,nhess;
47 double hmax;
48 char hatext[60];
49
50 int *keep; // keep[MAXATOM];
51
52 keep = ivector(0,natom);
53
54 nhess = 0;
55 for (i=1; i <= natom; i++)
56 {
57
58 hinit[nhess] = 0;
59 hstop[nhess] = 0;
60 hdiag[nhess] = 0.0;
61 nhess++;
62 hinit[nhess] = 0;
63 hstop[nhess] = 0;
64 hdiag[nhess] = 0.0;
65 nhess++;
66 hinit[nhess] = 0;
67 hstop[nhess] = 0;
68 hdiag[nhess] = 0.0;
69 nhess++;
70 }
71 /* periodic boundary conditions set here */
72
73 // if (pot.use_picalc) piseq(0,0);
74
75 nhess= 0;
76 for (i=1; i <= natom; i++)
77 {
78 for(k=1; k <= natom; k++)
79 {
80 for(j=0; j < 3; j++)
81 {
82 hess.hessx[k][j] = 0.0;
83 hess.hessy[k][j] = 0.0;
84 hess.hessz[k][j] = 0.0;
85 }
86 }
87
88 if (atom[i].use)
89 {
90 if (pot.use_bond)ebond2(i);
91 if (pot.use_angle)eangle2(i);
92 if (pot.use_strbnd)estrbnd2(i);
93 if (pot.use_opbend_wilson) eopbend_wilson2(i);
94 if (pot.use_tors)etorsion2(i);
95
96 if (pot.use_hal)ehal2(i);
97 if (pot.use_bufcharge)ebufcharge2(i);
98
99 if (pot.use_geom) egeom2(i);
100
101 hdiag[(i-1)*3] += hess.hessx[i][0];
102 hdiag[(i-1)*3+1] += hess.hessy[i][1];
103 hdiag[(i-1)*3+2] += hess.hessz[i][2];
104
105 for (k=i+1; k <= natom; k++)
106 {
107 keep[k] = FALSE;
108 if (atom[k].use)
109 {
110 hmax = fabs(hess.hessx[k][0]);
111 for (j=0; j < 3;j++)
112 {
113 if (fabs(hess.hessx[k][j]) > hmax)
114 hmax = fabs(hess.hessx[k][j]);
115 }
116 for (j=0; j < 3;j++)
117 {
118 if (fabs(hess.hessy[k][j]) > hmax)
119 hmax = fabs(hess.hessy[k][j]);
120 }
121 for (j=0; j < 3;j++)
122 {
123 if (fabs(hess.hessz[k][j]) > hmax)
124 hmax = fabs(hess.hessz[k][j]);
125 }
126 if (hmax >= hesscut) keep[k] = TRUE;
127 }
128 }
129
130 hinit[(i-1)*3] = nhess; // x component of new atom
131
132 hindex[nhess] = 3*(i-1) + 1;
133 h[nhess] = hess.hessx[i][1];
134 nhess++;
135 hindex[nhess] = 3*(i-1) + 2;
136 h[nhess] = hess.hessx[i][2];
137 nhess++;
138
139 for (k=i+1; k <= natom; k++)
140 {
141 if (keep[k])
142 {
143 for (j=0; j < 3; j++)
144 {
145 hindex[nhess] = 3*(k-1) +j;
146 h[nhess] = hess.hessx[k][j];
147 nhess++;
148 }
149 }
150 }
151 hstop[(i-1)*3] = nhess;
152
153 hinit[(i-1)*3+1] = nhess; // y component of atom i
154 hindex[nhess] = 3*(i-1) + 2;
155 h[nhess] = hess.hessy[i][2];
156 nhess++;
157 for (k=i+1; k <= natom; k++)
158 {
159 if (keep[k])
160 {
161 for (j=0; j < 3; j++)
162 {
163 hindex[nhess] = 3*(k-1) +j;
164 h[nhess] = hess.hessy[k][j];
165 nhess++;
166 }
167 }
168 }
169 hstop[(i-1)*3+1] = nhess;
170 hinit[(i-1)*3+2] = nhess;
171 for (k=i+1; k <= natom; k++)
172 {
173 if (keep[k])
174 {
175 for (j=0; j < 3; j++)
176 {
177 hindex[nhess] = 3*(k-1) +j;
178 h[nhess] = hess.hessz[k][j];
179 nhess++;
180 }
181 }
182 }
183 hstop[(i-1)*3+2] = nhess;
184
185 if (nhess > maxhess)
186 {
187 sprintf(hatext,"Too many hessian elements: %d elements %d max\n",nhess,maxhess);
188 message_alert(hatext,"Error");
189 fprintf(pcmlogfile,"Too many hessian elements !\n");
190 strcpy(Openbox.fname,"Hesserr.pcm");
191 pcmfout(2);
192 exit(0);
193 }
194 }
195 }
196 free_ivector(keep, 0,natom);
197 }