ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin2.c
Revision: 58
Committed: Mon Dec 1 06:44:59 2008 UTC (13 years, 5 months ago) by wdelano
Original Path: trunk/src/mengine/src/pcmwin2.c
File size: 38740 byte(s)
Log Message:
code merge 20081130
Line File contents
1 #define EXTERN extern
2 #include "pcwin.h"
3 #include "pcmod.h"
4 #include "field.h"
5 #include "atom_k.h"
6
7 void deleteatom(int);
8 int make_atom(int, float, float, float,char *);
9 void make_bond(int, int, int);
10 void hadd(void);
11 void hdel(int);
12 void type(void);
13 void set_atomtypes(int);
14 void nhadd(void);
15 void hcoord(void);
16 void revec(void);
17 void xyplan(float *, float *, float *);
18 void xaxis(float *,float *,float *);
19 void getvec(void);
20 void reseq(void);
21 void generate_bonds(void);
22 void set_atomdata(int, int,int,int,int,int);
23 void message_alert(char *, char *);
24
25 EXTERN struct t_minim_control {
26 int type, method, field, added_const;
27 char added_path[256],added_name[256];
28 } minim_control;
29
30 void hadd(void)
31 {
32 type();
33 if (field.type == MM3)
34 {
35 nhadd();
36 hcoord();
37 hdel(1);
38 } else if (field.type == AMBER)
39 {
40 nhadd();
41 hcoord();
42 hdel(1);
43 } else if (field.type == OPLSAA)
44 {
45 nhadd();
46 hcoord();
47 hdel(1);
48 } else if (field.type == MMFF94)
49 {
50 nhadd();
51 hcoord();
52 hdel(1);
53 } else
54 {
55 set_atomtypes(MMX);
56 nhadd();
57 hcoord();
58 }
59 type();
60 generate_bonds();
61
62 }
63
64 void nhadd()
65 {
66 int i, i1, ik, im, it, itads, j, jk, nhadds, newatom;
67 int jji, bo,ncount;
68 static int iadd[MAXVDW]={ 0,0,0,0,0,3,1,6,0,1,
69 0,0,0,0,3,6,3,0,0,0,
70 0,0,0,0,6,0,0,0,0,0,
71 0,0,0,3,0,0,6,1,1,0,
72 0,4,0,0,0,6,0,6,0,0,
73 3,3,3,0,6,0,0,0,0,0,
74 0,0,0,0,0,1,6,0,0,0,
75 0,0,0,0,0};
76
77 /* this routine calculates the #'s of h to be added
78 * to carbon and packs it symbolically into iat(i,j) table
79 * as well as into the connection table via the appropriate
80 * manipulation instructions. hydrogen coordinates
81 * are not added by this routine. lp are also added to oxygen. */
82 /* 0 to disregard : 1 and 4 for lp's only
83 * 2 and 5 for h's only and 3 and 6 for h and lp
84 * add 3 for trigonal atoms */
85
86 ncount = natom;
87 for( i = 1 ; i <= ncount; i++ )
88 {
89 if( atom[i].mmx_type > 0 && (atom[i].mmx_type < 300) && atom[i].atomnum != 1 )
90 {
91 jji = 0;
92 bo = 0;
93 for (j=0; j < MAXIAT; j++)
94 {
95 if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
96 {
97 jji ++;
98 bo += atom[i].bo[j];
99 }
100 }
101 nhadds = atom_def.ligands[atom[i].mmx_type] - jji;
102
103 if (atom[i].mmx_type >= 11 && atom[i].mmx_type <= 14)
104 nhadds = 0;
105 if (atom[i].mmx_type == 41 || atom[i].mmx_type == 46)
106 nhadds = atom_def.ligands[atom[i].mmx_type] - bo;
107 if (atom[i].mmx_type == 49 || atom[i].mmx_type == 50 || atom[i].mmx_type == 51)
108 nhadds = atom_def.ligands[atom[i].mmx_type] - bo;
109
110 /* add atom and bond it to atom i */
111 if( nhadds > 0 )
112 {
113 for( im = 1; im <= nhadds; im++ )
114 {
115 newatom = make_atom(5,0.F,0.F,0.F,"H");
116 set_atomdata(newatom,5,5,5,34,29);
117 if (newatom == -1)
118 {
119 message_alert("Error in nhadd newatom = -1","Error");
120 hdel(0);
121 return;
122 }
123 make_bond(i,newatom,1);
124 }
125 }
126 }
127 }
128
129 /* *** check for lp's *** */
130 if (field.type == MMX)
131 {
132 ncount = natom;
133 for( i = 1 ; i <= ncount; i++ )
134 {
135 it = iadd[atom[i].mmx_type-1];
136 if( it >= 1 )
137 {
138 itads = atom[i].mmx_type ;
139 if (it == 1 || it == 3 || it == 4 || it == 6 )
140 {
141 nhadds = 4;
142 for( i1 = 0; i1 < MAXIAT; i1++ )
143 {
144 if( atom[i].bo[i1] != 0 && atom[i].bo[i1] != 9 )
145 nhadds = nhadds - atom[i].bo[i1];
146 }
147 if( itads == 6 || itads == 15 || itads == 34 || itads == 42 || itads == 53 || itads == 48 )
148 {
149 for( ik = 0; ik < MAXIAT; ik++ )
150 {
151 if( atom[i].iat[ik] != 0 && atom[i].bo[ik] != 9 )
152 {
153 jk = atom[atom[i].iat[ik]].mmx_type;
154 if( jk == 2 || jk == 3 || jk == 4 || jk == 9
155 || jk == 10 || jk == 30 || jk == 37 || jk == 40 )
156 {
157 nhadds = nhadds - 1;
158 break;
159 }
160 }
161 }
162 }
163
164 /* correction for sulfoxide */
165 if( itads == 17 )
166 nhadds = nhadds + 1;
167
168 if (itads == 66)
169 nhadds = 2;
170
171 if( nhadds > 0 )
172 {
173 for( i1 = 1; i1 <= nhadds; i1++ )
174 {
175 newatom = make_atom(20, 0.F,0.F,0.F,"");
176 set_atomdata(newatom,20,0,0,0,0);
177 make_bond(i,newatom,1);
178 }
179 }
180 }
181 }
182 }
183 }
184 return;
185 }
186 // ==================================
187 void hcoord()
188 {
189 int i, i1, i3, i4, i5, i6, i7, i8,
190 i9, iadj[10], ii, ii1, ii2, ii3, iii,
191 ij, ijk[3], it, j, j1, j2, j3, j4, j5, j8,
192 j9, k1, k2, l1, na, nhadds, ni;
193 float adj, dis, xpnt, ypnt, zpnt;
194
195 getvec();
196 for( i = 1; i <= natom; i++ )
197 {
198 it = atom[i].mmx_type;
199 if( it < 300)
200 {
201 for( j = 0; j < 10; j++ )
202 iadj[j] = 0;
203
204 for( j = 0; j < 3; j++ )
205 ijk[j] = 0;
206
207 nhadds = 0;
208 l1 = 0;
209 for( j = 0; j < MAXIAT; j++ )
210 {
211 if( atom[i].iat[j] != 0 && atom[i].bo[j] != 9 )
212 {
213 if( atom[atom[i].iat[j]].mmx_type == 5 || atom[atom[i].iat[j]].mmx_type == 20 )
214 nhadds++ ;
215 else
216 {
217 iadj[l1] = atom[i].iat[j];
218 l1++;
219 }
220 }
221 }
222 if( nhadds != 0 )
223 {
224 /* **** calculate the coordinates of the center we
225 * are going to move to the x-axis */
226 xpnt = atom[i].x;
227 ypnt = atom[i].y;
228 zpnt = atom[i].z;
229 for( ii = 1; ii <= (natom + 4); ii++ )
230 {
231 atom[ii].x = atom[ii].x - xpnt;
232 atom[ii].y = atom[ii].y - ypnt;
233 atom[ii].z = atom[ii].z - zpnt;
234 }
235 xpnt = 0.F;
236 ypnt = 0.F;
237 zpnt = 0.F;
238 adj = 0.F;
239 for( i1 = 0; i1 < 10; i1++ )
240 {
241 if( iadj[i1] != 0 )
242 {
243 na = iadj[i1];
244 adj += 1.0001F ;
245 xpnt += atom[na].x;
246 ypnt += atom[na].y;
247 zpnt += atom[na].z;
248 }
249 }
250 if( adj > 0.01F )
251 {
252 xpnt /= adj;
253 ypnt /= adj;
254 zpnt /= adj;
255 xaxis( &xpnt, &ypnt, &zpnt );
256 if( xpnt > 0.F )
257 {
258 for( i3 = 1; i3 <= (natom + 4); i3++ )
259 atom[i3].x = -atom[i3].x;
260 }
261 /* **** after moving, check the # of h added to the carbon i */
262 if( nhadds == 1 )
263 {
264 /* **** there is only 1 h added to the carbon i */
265 for( i4 = 0; i4 < MAXIAT; i4++ )
266 {
267 if( atom[i].iat[i4] != 0 && atom[i].bo[i4] != 9 )
268 {
269 i5 = atom[i].iat[i4];
270 if( atom[i5].mmx_type == 5 || atom[i5].mmx_type == 20 )
271 {
272 /* i5 is the atom to attach */
273 if (adj < 2.0) // fix for OH in MMFF
274 {
275 atom[i5].x = 0.70;
276 atom[i5].y = 0.70;
277 atom[i5].z = 0.00;
278 } else
279 {
280 atom[i5].x = 1.10F;
281 if( atom[i5].mmx_type == 20 )
282 atom[i5].x = 0.6F;
283 atom[i5].y = 0.00F;
284 atom[i5].z = 0.00F;
285 }
286 if( atom[i].atomnum == 16 )
287 {
288 atom[i5].x = 0.F;
289 atom[i5].y = 0.3796F;
290 atom[i5].z = 1.043F;
291 }
292 }
293 }
294 }
295 }
296 /* **** there are 2 h added to the carbon i */
297 if( nhadds == 2 )
298 {
299 if( it == 2 || it == 3 || it == 29 || it == 30 || it == 9 || it == 41 || it == 46
300 || it == 48 || it == 7 || it == 37 || it == 39
301 || it == 42 || it == 38 || it == 40 || it == 44 || it == 66)
302 {
303 for( i8 = 0; i8 < 10; i8++ )
304 {
305 if( iadj[i8] != 0 )
306 {
307 ii1 = iadj[i8];
308 for( ii2 = 0; ii2 < MAXIAT; ii2++ )
309 {
310 ii3 = 0;
311 if( atom[ii1].iat[ii2] != 0 && atom[ii1].iat[ii2] != i && atom[ii1].bo[ii2] != 9 )
312 {
313 ii3 = atom[ii1].iat[ii2];
314 xpnt = atom[ii3].x;
315 ypnt = atom[ii3].y;
316 zpnt = atom[ii3].z;
317 xyplan( &xpnt, &ypnt, &zpnt );
318 ij = 0;
319 for( i9 = 0; i9 < MAXIAT; i9++ )
320 {
321 if( atom[i].iat[i9] != 0 && atom[i].bo[i9] != 9 )
322 {
323 iii = atom[i].iat[i9];
324 if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 )
325 {
326 ijk[ij] = iii;
327 ij++;
328 if( ij == 2 )
329 {
330 /* *** =ch2 *** */
331 atom[ijk[0]].x = 0.5582F;
332 atom[ijk[1]].x = 0.5582F;
333 atom[ijk[0]].y = 0.9478F;
334 atom[ijk[1]].y = -0.9478F;
335 atom[ijk[0]].z = 0.00F;
336 atom[ijk[1]].z = 0.00F;
337 if( it == 7 || it == 38 || it == 42 || it == 39 || it == 66)
338 {
339 atom[ijk[0]].x = 0.30F;
340 atom[ijk[1]].x = 0.30F;
341 atom[ijk[0]].y = 0.52F;
342 atom[ijk[1]].y = -0.52F;
343 }
344 if( atom[i].mmx_type == 37 )
345 {
346 atom[ijk[1]].x = 0.35F;
347 atom[ijk[1]].y = -0.5F;
348 }
349 }
350 }
351 }
352 }
353 break;
354 }
355 if (ii3 == 0) // C=O with no lone pairs
356 {
357 ij = 0;
358 for( i9 = 0; i9 < MAXIAT; i9++ )
359 {
360 if( atom[i].iat[i9] != 0 && atom[i].bo[i9] != 9 )
361 {
362 iii = atom[i].iat[i9];
363 if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 )
364 {
365 ijk[ij] = iii;
366 ij++;
367 if( ij == 2 )
368 {
369 /* *** =ch2 *** */
370 atom[ijk[0]].x = 0.5582F;
371 atom[ijk[1]].x = 0.5582F;
372 atom[ijk[0]].y = 0.9478F;
373 atom[ijk[1]].y = -0.9478F;
374 atom[ijk[0]].z = 0.00F;
375 atom[ijk[1]].z = 0.00F;
376 if( it == 7 || it == 38 || it == 42 || it == 39 || it == 66)
377 {
378 atom[ijk[0]].x = 0.30F;
379 atom[ijk[1]].x = 0.30F;
380 atom[ijk[0]].y = 0.52F;
381 atom[ijk[1]].y = -0.52F;
382 }
383 if( atom[i].mmx_type == 37 )
384 {
385 atom[ijk[1]].x = 0.35F;
386 atom[ijk[1]].y = -0.5F;
387 }
388 }
389 }
390 }
391 }
392 break;
393 }
394 }
395 }
396 }
397 } else
398 {
399 /* ** there are 2 atoms attached to carbon i,
400 * one of which is "ni" and move these two atoms to xy-plane */
401 for( i6 = 0; i6 < 10; i6++ )
402 {
403 if( iadj[i6] != 0 )
404 {
405 ni = iadj[i6];
406 xpnt = atom[ni].x;
407 ypnt = atom[ni].y;
408 zpnt = atom[ni].z;
409 xyplan( &xpnt, &ypnt, &zpnt );
410 }
411 }
412 ij = 0;
413 for( i7 = 0; i7 < MAXIAT; i7++ )
414 {
415 if( atom[i].iat[i7] != 0 && atom[i].bo[i7] != 9 )
416 {
417 iii = atom[i].iat[i7];
418 if( atom[iii].mmx_type == 5 || atom[iii].mmx_type == 20 )
419 {
420 ijk[ij] = iii;
421 ij++ ;
422 if( ij == 2 )
423 {
424 /* *** -ch2- type 1 atom *** */
425 atom[ijk[0]].x = 0.6758F;
426 atom[ijk[1]].x = 0.6758F;
427 atom[ijk[0]].y = 0.00F;
428 atom[ijk[1]].y = 0.00F;
429 atom[ijk[0]].z = 0.8807F;
430 atom[ijk[1]].z = -0.8807F;
431 if( atom[i].mmx_type == 25 || atom[i].mmx_type == 46 || atom[i].mmx_type == 67)
432 {
433 atom[ijk[0]].x = 0.3796F;
434 atom[ijk[1]].x = 0.3796F;
435 atom[ijk[0]].y = 1.043F;
436 atom[ijk[1]].y = -0.5215F;
437 atom[ijk[0]].z = 0.0000F;
438 atom[ijk[1]].z = 0.9032F;
439 }
440 if( atom[i].mmx_type == 22 || atom[i].mmx_type == 30 )
441 {
442 /* *** -ch2- type 22 and 3o atoms *** */
443 atom[ijk[0]].x = 0.2F;
444 atom[ijk[1]].x = 0.2F;
445 atom[ijk[0]].y = 0.0F;
446 atom[ijk[1]].y = 0.0F;
447 atom[ijk[0]].z = 1.0F;
448 atom[ijk[1]].z = -1.0F;
449 }
450 /* *** -o(lp)2- *** */
451 if (atom[i].mmx_type == 6 && (atom[iadj[0]].mmx_type == 3 ||
452 atom[iadj[0]].mmx_type == 2 || atom[iadj[0]].mmx_type == 30 ||
453 atom[iadj[0]].mmx_type == 37 || atom[iadj[0]].mmx_type == 40) )
454 {
455 atom[ijk[0]].x = 0.450F;
456 atom[ijk[1]].x = 0.2488F;
457 atom[ijk[0]].z = 0.0F;
458 atom[ijk[1]].z = 0.0F;
459 atom[ijk[0]].y = 0.9500F;
460 atom[ijk[1]].y = -0.5460F;
461 }
462 else if( atom[i].mmx_type == 6 || atom[i].mmx_type == 15 || atom[i].mmx_type ==
463 42 || atom[i].mmx_type == 34 || atom[i].mmx_type == 35 || atom[i].mmx_type == 53 )
464 {
465 atom[ijk[0]].x = 0.2488F;
466 atom[ijk[1]].x = 0.2488F;
467 atom[ijk[0]].y = 0.0F;
468 atom[ijk[1]].y = 0.0F;
469 atom[ijk[0]].z = 0.5460F;
470 atom[ijk[1]].z = -0.5460F;
471 }
472 }
473 }
474 }
475 }
476 }
477 }
478 /* end of nhadds = 2
479 * **** there are 3 h added to the carbon i *** */
480 if( nhadds == 3 )
481 {
482 for( j1 = 0; j1 < 10; j1++ )
483 {
484 if( iadj[j1] != 0 )
485 {
486 j2 = iadj[j1];
487 if( atom[j2].atomnum != 1 && atom[j2].mmx_type != 20 )
488 na = j2;
489 if( atom[j2].mmx_type == 2 || atom[j2].mmx_type == 3 )
490 {
491 for( j3 = 0; j3 < MAXIAT; j3++ )
492 {
493 if( (atom[j2].iat[j3] != 0 && atom[j2].iat[j3] != i) &&
494 atom[j2].bo[j3] != 9 )
495 {
496 j4 = atom[j2].iat[j3];
497 if( atom[j4].mmx_type == 2 || (atom[j4].mmx_type == 7 && atom[j4].mmx_type == 38) )
498 {
499 xpnt = atom[j4].x;
500 ypnt = atom[j4].y;
501 zpnt = atom[j4].z;
502 xyplan( &xpnt, &ypnt, &zpnt );
503 if( atom[j4].y < 0.F )
504 {
505 for( j5 = 1; j5 <= (natom + 4); j5++ )
506 atom[j5].y = -atom[j5].y;
507 }
508 }
509 }
510 }
511 }else
512 {
513 /* end of type 2 and type 3 adjacent* */
514 for( j8 = 0; j8 < MAXIAT; j8++ )
515 {
516 if( (atom[j2].iat[j8] != 0 && atom[j2].iat[j8] != i)
517 && atom[j2].bo[j8] != 9 )
518 {
519 j9 = atom[j2].iat[j8];
520 // if( atom[j9].mmx_type != 5 && atom[j9].mmx_type != 20 )
521 if(atom[j9].mmx_type != 20 )
522 {
523 xpnt = atom[j9].x;
524 ypnt = atom[j9].y;
525 zpnt = atom[j9].z;
526 xyplan( &xpnt, &ypnt, &zpnt );
527 if( atom[j9].y > 0.F )
528 {
529 for( k1 = 1; k1 <= (natom + 4); k1++ )
530 atom[k1].y = -atom[k1].y;
531 }
532 }
533 }
534 }
535 }
536 /* do type 4 here */
537 ij = 0;
538 for( k2 = 0; k2 < MAXIAT; k2++ )
539 {
540 if( !(atom[i].iat[k2] == 0 || atom[i].bo[k2] == 9) )
541 {
542 iii = atom[i].iat[k2];
543 if( !(atom[iii].mmx_type != 5 && atom[iii].mmx_type != 20) )
544 {
545 ijk[ij] = iii;
546 ij++;
547 if( ij == 3 )
548 {
549 dis = 0.F;
550 /* *** -ch3 ***
551 * *** add the hydrogens displaced along x enough
552 * to give a normal length c-c bond *** */
553 atom[ijk[0]].x = 0.3796F + dis;
554 atom[ijk[1]].x = 0.3796F + dis;
555 atom[ijk[2]].x = 0.3796F + dis;
556 atom[ijk[0]].y = 1.043F;
557 atom[ijk[1]].y = -0.5215F;
558 atom[ijk[2]].y = -0.5215F;
559 atom[ijk[0]].z = 0.00F;
560 atom[ijk[1]].z = 0.9032F;
561 atom[ijk[2]].z = -0.9032F;
562 }
563 }
564 }
565 }
566 }
567 }
568 }
569 }
570 }
571 }
572 }
573 /* *** reorient the structure in starting orientation and
574 * turn off atom markers *** */
575 revec();
576 return;
577 }
578 // ========================================================
579 void revec()
580 {
581 int i;
582 float xpnt, ypnt, zpnt;
583 /* *** this routine used the unit vectors and origin saved
584 * by 'subroutine getvec' to return the molecule to
585 * its original orientation *** */
586
587 if( natom <= MAXATOM - 4 )
588 {
589 for( i = 1; i <= (natom + 4); i++ )
590 {
591 atom[i].x = atom[i].x - atom[natom + 4].x;
592 atom[i].y = atom[i].y - atom[natom + 4].y;
593 atom[i].z = atom[i].z - atom[natom + 4].z;
594 }
595 xpnt = atom[natom + 1].x;
596 ypnt = atom[natom + 1].y;
597 zpnt = atom[natom + 1].z;
598 xaxis( &xpnt, &ypnt, &zpnt );
599 xpnt = atom[natom + 2].x;
600 ypnt = atom[natom + 2].y;
601 zpnt = atom[natom + 2].z;
602 xyplan( &xpnt, &ypnt, &zpnt );
603 if( atom[natom + 1].x < 0. )
604 {
605 for( i = 1; i <= (natom + 4); i++ )
606 atom[i].x = -atom[i].x;
607 }
608 if( atom[natom + 2].y < 0. )
609 {
610 for( i = 1; i <= (natom + 4); i++ )
611 atom[i].y = -atom[i].y;
612 }
613 if( atom[natom + 3].z < 0. )
614 {
615 for( i = 1; i <= (natom + 4); i++ )
616 atom[i].z = -atom[i].z;
617 }
618 }
619 return;
620 }
621 /* ------------------------ */
622 void xyplan(float *xr, float *yr, float *zr)
623 {
624 int i;
625 float cos3, denom, savex, sin3;
626
627 /* this subroutine rotates the molecule to place atom (xr,yr,*zr)
628 * in the x,y-plane. note0 y-coordinate will be +. */
629
630 denom = sqrt( *yr**yr + *zr**zr );
631 if( fabs( denom ) >= .00001F )
632 {
633 sin3 = *zr/denom;
634 cos3 = *yr/denom;
635 for( i = 1; i <= (natom + 4); i++ )
636 {
637 savex = atom[i].y;
638 atom[i].y = atom[i].y*cos3 + atom[i].z*sin3;
639 atom[i].z = atom[i].z*cos3 - savex*sin3;
640 }
641 savex = *yr;
642 *yr = *yr*cos3 + *zr*sin3;
643 *zr = *zr*cos3 - savex*sin3;
644 }
645 return;
646 }
647 /* ------------------------ */
648 void xaxis(float *xpnt,float *ypnt,float *zpnt)
649 {
650 int i;
651 float cos1, cos2, denom, save, sin1, sin2;
652
653 /* this subroutine rotates the molecule about z-axis and then y-axis to
654 * place point with coordinates (i*xpnt,iypnt,izpnt) along x-axis. */
655
656 denom = sqrt( *xpnt**xpnt + *ypnt**ypnt );
657 if( fabs( denom ) < .00001F )
658 {
659 sin1 = 0.F;
660 cos1 = 1.F;
661 }else
662 {
663 sin1 = *ypnt/denom;
664 cos1 = *xpnt/denom;
665 for( i = 1; i <= (natom + 4); i++ )
666 {
667 save = atom[i].x;
668 atom[i].x = atom[i].x*cos1 + atom[i].y*sin1;
669 atom[i].y = atom[i].y*cos1 - save*sin1;
670 }
671 save = *xpnt;
672 *xpnt = *xpnt*cos1 + *ypnt*sin1;
673 *ypnt = *ypnt*cos1 - save*sin1;
674 }
675 denom = sqrt( *xpnt**xpnt + *zpnt**zpnt );
676 if( fabs( denom ) >= .00001F )
677 {
678 sin2 = *zpnt/denom;
679 cos2 = *xpnt/denom;
680 for( i = 1; i <= (natom + 4); i++ )
681 {
682 save = atom[i].x;
683 atom[i].x = atom[i].x*cos2 + atom[i].z*sin2;
684 atom[i].z = atom[i].z*cos2 - save*sin2;
685 }
686 save = *xpnt;
687 *xpnt = *xpnt*cos2 + *zpnt*sin2;
688 *zpnt = *zpnt*cos2 - save*sin2;
689 }
690 return;
691 }
692 /* ============================================== */
693 void getvec()
694 {
695
696 /* *** this routine saves unit vectors and origin in the
697 * n+1 - n+4 slots of x, y, and z *** */
698
699 if(natom > MAXATOM - 4 )
700 return;
701 atom[natom + 1].x = 1.F;
702 atom[natom + 1].y = 0.F;
703 atom[natom + 1].z = 0.F;
704 atom[natom + 2].x = 0.F;
705 atom[natom + 2].y = 1.F;
706 atom[natom + 2].z = 0.F;
707 atom[natom + 3].x = 0.F;
708 atom[natom + 3].y = 0.F;
709 atom[natom + 3].z = 1.F;
710 atom[natom + 4].x = 0.F;
711 atom[natom + 4].y = 0.F;
712 atom[natom + 4].z = 0.F;
713 return;
714 }
715 /* ============================================== */
716 void hdel(int lptest)
717 {
718 /* lptest = 0 to remove H and lone pair */
719 /* lptest = 1 to remove lone pairs only */
720 /* lptest = 2 to remove lp & reset various groups for amber */
721 int i, iatt, iatt1, iatyp, it, it1,ntemp;
722 int j, jatom;
723 /* *** delete hydrogens *** */
724 ntemp = natom;
725 for( i = 1; i <= ntemp; i++ )
726 {
727 if( atom[i].mmx_type != 0 )
728 {
729 iatyp = atom[i].mmx_type;
730 if( (lptest == 0 && (iatyp == 5 || iatyp == 20 )) || (lptest >= 1 && iatyp == 20) )
731 {
732 iatt = atom[i].iat[0];
733 iatt1 = atom[i].iat[1];
734 if( iatt != 0 || iatt1 != 0)
735 {
736 it = atom[iatt].mmx_type;
737 it1 = atom[iatt1].mmx_type;
738 if (it >= 11 && it <= 14) // halogens
739 goto L_10;
740 if( (it < 300) && (it1 < 300)) // attached metals
741 deleteatom(i);
742 }
743 }
744 }
745 L_10:
746 continue;
747 }
748 reseq();
749 generate_bonds();
750 return;
751 }
752 /* -------------------------- */
753 void reseq()
754 {
755 int i, ia, ib, iplus, ixx, j, jincr, k;
756
757 /* ***this subroutine will repack the atom
758 * ***connection table, making sure to fill in
759 * **all "holes" caused by deleting atoms or bonds.
760 * ***the necessary common blocks involve arrays
761 * of the connection table (contb) */
762
763 for( i = 1; i < MAXATOM; i++ )
764 {
765 if( atom[i].type != 0 )
766 {
767 for( ia = 0; ia < (MAXIAT - 1); ia++ )
768 {
769 if( atom[i].iat[ia] == 0 )
770 {
771 /* *** for every zero entry in the array iat
772 * for atom i, want to move down all remaining
773 * entries of iat by one slot; after doing so
774 * want to store that zero entry in the last
775 * slot--iat *** */
776 for( ib = ia + 1; ib < MAXIAT; ib++ )
777 {
778 atom[i].iat[ib - 1] = atom[i].iat[ib];
779 atom[i].bo[ib - 1] = atom[i].bo[ib];
780 }
781 atom[i].iat[MAXIAT - 1] = 0;
782 }
783 }
784 }
785 }
786
787 /* *** now want to renumber atoms to leave
788 * no "holes" caused by the deleted atoms *** */
789 jincr = 0;
790 /* i points to first zero entry
791 * kincr points to next real entry */
792 for( i = 1; i < (MAXATOM - 1); i++ )
793 {
794 iplus = 0;
795 if( atom[i].type == 0 )
796 {
797 for( jincr = i; jincr < MAXATOM; jincr++ )
798 {
799 if( atom[jincr].type != 0 )
800 {
801 iplus = jincr;
802 break;
803 }
804 }
805 if (iplus == 0)
806 break;
807 for( j = 0; j < MAXIAT; j++ )
808 {
809 atom[i].iat[j] = atom[iplus].iat[j];
810 atom[i].bo[j] = atom[iplus].bo[j];
811 }
812
813 atom[i].x = atom[iplus].x;
814 atom[i].y = atom[iplus].y;
815 atom[i].z = atom[iplus].z;
816 atom[i].type = atom[iplus].type;
817 atom[i].tclass = atom[iplus].tclass;
818 atom[i].type = atom[iplus].type;
819 atom[i].mmx_type = atom[iplus].mmx_type;
820 atom[i].mm3_type = atom[iplus].mm3_type;
821 atom[i].amber_type = atom[iplus].amber_type;
822 atom[i].mmff_type = atom[iplus].mmff_type;
823 atom[i].atomnum = atom[iplus].atomnum;
824 atom[i].formal_charge = atom[iplus].formal_charge;
825 atom[i].atomwt = atom[iplus].atomwt;
826 atom[i].use = atom[iplus].use;
827 strcpy(atom[i].name,atom[iplus].name);
828 atom[i].charge = atom[iplus].charge;
829 atom[i].flags = atom[iplus].flags;
830 atom[i].radius = atom[iplus].radius;
831 atom[i].energy = atom[iplus].energy;
832 /* *** delete the atom whose statistics were just
833 * stored in the i-th atom's connection table *** */
834 for( ixx = 0; ixx < MAXIAT; ixx++ )
835 {
836 atom[iplus].iat[ixx] = 0;
837 atom[iplus].bo[ixx] = 0;
838 }
839 atom[iplus].x = 0.0F;
840 atom[iplus].y = 0.0F;
841 atom[iplus].z = 0.0F;
842 atom[iplus].type = 0;
843 atom[iplus].tclass = 0;
844 atom[iplus].mmx_type = 0;
845 atom[iplus].mm3_type = 0;
846 atom[iplus].amber_type = 0;
847 atom[iplus].mmff_type = 0;
848 atom[iplus].atomnum = 0;
849 atom[iplus].formal_charge = 0;
850 atom[iplus].atomwt = 0.0F;
851 atom[iplus].use = 0;
852 atom[iplus].charge = 0.0F;
853 atom[iplus].flags = 0;
854 atom[iplus].radius = 0.0F;
855 atom[iplus].energy = 0.0F;
856 strcpy(atom[iplus].name,"");
857
858 for( j = 0; j < MAXIAT; j++ )
859 {
860 if( atom[i].iat[j] != 0 )
861 {
862 for( k = 0; k < MAXIAT; k++ )
863 {
864 if( atom[atom[i].iat[j]].iat[k] == iplus )
865 atom[atom[i].iat[j]].iat[k] = i;
866 }
867 }
868 }
869 }
870 }
871 /* ***now want to do the same resequencing for iat array */
872
873 /* recalc number of atoms */
874 natom = 0;
875 for (i = 1; i < MAXATOM; i++)
876 {
877 if (atom[i].type != 0)
878 natom++;
879 }
880 return;
881 }