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

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