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

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 message_alert(char *,char *);
42
43 void hessian(int maxhess, double *h, int *hinit, int *hstop,int *hindex, double *hdiag)
44 {
45 int i,j,k,nhess;
46 double hmax;
47 char hatext[60];
48
49 int *keep; // keep[MAXATOM];
50
51 keep = ivector(0,natom);
52
53 nhess = 0;
54 for (i=1; i <= natom; i++)
55 {
56
57 hinit[nhess] = 0;
58 hstop[nhess] = 0;
59 hdiag[nhess] = 0.0;
60 nhess++;
61 hinit[nhess] = 0;
62 hstop[nhess] = 0;
63 hdiag[nhess] = 0.0;
64 nhess++;
65 hinit[nhess] = 0;
66 hstop[nhess] = 0;
67 hdiag[nhess] = 0.0;
68 nhess++;
69 }
70 /* periodic boundary conditions set here */
71
72 // if (pot.use_picalc) piseq(0,0);
73
74 nhess= 0;
75 for (i=1; i <= natom; i++)
76 {
77 for(k=1; k <= natom; k++)
78 {
79 for(j=0; j < 3; j++)
80 {
81 hess.hessx[k][j] = 0.0;
82 hess.hessy[k][j] = 0.0;
83 hess.hessz[k][j] = 0.0;
84 }
85 }
86
87 if (atom[i].use)
88 {
89 if (pot.use_bond)ebond2(i);
90 if (pot.use_angle)eangle2(i);
91 if (pot.use_strbnd)estrbnd2(i);
92 if (pot.use_opbend_wilson) eopbend_wilson2(i);
93 if (pot.use_tors)etorsion2(i);
94
95 if (pot.use_hal)ehal2(i);
96 if (pot.use_bufcharge)ebufcharge2(i);
97
98 hdiag[(i-1)*3] += hess.hessx[i][0];
99 hdiag[(i-1)*3+1] += hess.hessy[i][1];
100 hdiag[(i-1)*3+2] += hess.hessz[i][2];
101
102 for (k=i+1; k <= natom; k++)
103 {
104 keep[k] = FALSE;
105 if (atom[k].use)
106 {
107 hmax = fabs(hess.hessx[k][0]);
108 for (j=0; j < 3;j++)
109 {
110 if (fabs(hess.hessx[k][j]) > hmax)
111 hmax = fabs(hess.hessx[k][j]);
112 }
113 for (j=0; j < 3;j++)
114 {
115 if (fabs(hess.hessy[k][j]) > hmax)
116 hmax = fabs(hess.hessy[k][j]);
117 }
118 for (j=0; j < 3;j++)
119 {
120 if (fabs(hess.hessz[k][j]) > hmax)
121 hmax = fabs(hess.hessz[k][j]);
122 }
123 if (hmax >= hesscut) keep[k] = TRUE;
124 }
125 }
126
127 hinit[(i-1)*3] = nhess; // x component of new atom
128
129 hindex[nhess] = 3*(i-1) + 1;
130 h[nhess] = hess.hessx[i][1];
131 nhess++;
132 hindex[nhess] = 3*(i-1) + 2;
133 h[nhess] = hess.hessx[i][2];
134 nhess++;
135
136 for (k=i+1; k <= natom; k++)
137 {
138 if (keep[k])
139 {
140 for (j=0; j < 3; j++)
141 {
142 hindex[nhess] = 3*(k-1) +j;
143 h[nhess] = hess.hessx[k][j];
144 nhess++;
145 }
146 }
147 }
148 hstop[(i-1)*3] = nhess;
149
150 hinit[(i-1)*3+1] = nhess; // y component of atom i
151 hindex[nhess] = 3*(i-1) + 2;
152 h[nhess] = hess.hessy[i][2];
153 nhess++;
154 for (k=i+1; k <= natom; k++)
155 {
156 if (keep[k])
157 {
158 for (j=0; j < 3; j++)
159 {
160 hindex[nhess] = 3*(k-1) +j;
161 h[nhess] = hess.hessy[k][j];
162 nhess++;
163 }
164 }
165 }
166 hstop[(i-1)*3+1] = nhess;
167 hinit[(i-1)*3+2] = nhess;
168 for (k=i+1; k <= natom; k++)
169 {
170 if (keep[k])
171 {
172 for (j=0; j < 3; j++)
173 {
174 hindex[nhess] = 3*(k-1) +j;
175 h[nhess] = hess.hessz[k][j];
176 nhess++;
177 }
178 }
179 }
180 hstop[(i-1)*3+2] = nhess;
181
182 if (nhess > maxhess)
183 {
184 sprintf(hatext,"Too many hessian elements: %d elements %d max\n",nhess,maxhess);
185 message_alert(hatext,"Error");
186 printf("Too many hessian elements !\n");
187 strcpy(Openbox.fname,"Hesserr.pcm");
188 pcmfout(2);
189 exit(0);
190 }
191 }
192 }
193 free_ivector(keep, 0,natom);
194 }