ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/pcmwin1.c
(Generate patch)
# Line 4 | Line 4
4   #include "field.h"
5   #include "atom_k.h"
6  
7 + double dihdrl(int i1,int i2,int i3,int i4);
8   void type_mmx(void);
8 void type_mm3(void);
9   void hdel(int);
10 void quick_type(void);
10   void set_atomdata(int,int,int,int,int,int);
11 < void message_alert(char *, char *);
12 < void set_atomtypes(int);        
11 >
12 > void set_atomtypes(int);
13 > void type(void);      
14 > static void quick_type(void);
15  
16   void set_atomtypes(int newtype)
17   {
# Line 198 | Line 199
199         message_alert("No typing rules for this force field are implemented!","ERROR");  
200   }
201   /* ------------------------ */
201
202   double dihdrl(int i1,int i2,int i3,int i4)
203   {
204          float           a1, a2, ang, b1, b2, c1, c2, siign,
# Line 252 | Line 252
252          return (dihdrl_v);
253   }                               /* end of function */
254   /* --------------------------------- */
255 void avgleg()
256 {
257        int i, iatyp, ib1, ib2, j;
258        int hasSulf;
259        float an, xx, xyzdis, yy, zz;
260
261        /****  if there is no bond, average=1.0 *** */
262
263        an = 0.0F;
264        xyzdis = 0.0F;
265        hasSulf = FALSE;
266        for( i = 1; i <= natom; i++ )
267        {
268                if( atom[i].type && atom[i].type != 20 )
269                {
270                   for( j = 0; j < MAXIAT; j++ )
271                   {
272                      if( atom[i].iat[j] && atom[i].bo[j] != 9 )
273                      {
274        /*             count bonds only once - don't count bonds to hydrogen */
275                          if( atom[i].iat[j] >= i )
276                          {
277                              ib1 = i;
278                              iatyp = atom[i].type;
279                              if (atom[i].atomnum == 16) hasSulf = TRUE;
280                              if(  iatyp != 5 && iatyp != 21 && iatyp != 23 && iatyp != 24 && iatyp != 28 && iatyp != 36 && iatyp != 15 && iatyp != 38 )
281                              {
282                                   ib2 = atom[i].iat[j];
283                                   iatyp = atom[atom[i].iat[j]].type;
284                                   if( iatyp != 20 )
285                                   {
286                                      if( iatyp != 5 && iatyp != 21 && iatyp != 23 && iatyp != 24 && iatyp != 28 && iatyp != 36 && iatyp != 15 && iatyp != 38)
287                                      {
288                                            an += 1.0F;
289                                            xx = atom[ib1].x - atom[ib2].x;
290                                            yy = atom[ib1].y - atom[ib2].y;
291                                            zz = atom[ib1].z - atom[ib2].z;
292                                            xyzdis += sqrt( xx*xx + yy*yy + zz*zz );
293                                       }
294                                    }
295                                }
296                           }
297                       }
298                   }
299                }
300        }
301        /* if an = 0 then there are no heavy atom bonds
302         *  need to redo calculations based on c-h bonds */
303        if( an != 0.0F )
304        {
305            if (hasSulf)
306                flags.avleg = 1.73F/(xyzdis/an);        /* assume normal c-c distance   */
307            else
308                flags.avleg = 1.53F/(xyzdis/an);        /* assume normal c-c distance   */
309        }
310        else if( natom > 1 )
311        {
312                xyzdis = 0.0F;
313                an = 0.0F;
314                for( i = 1; i <= natom; i++ )
315                {
316                        ib1 = i;
317                        for( j = 0; j < MAXIAT; j++ )
318                        {
319                                if( atom[i].iat[j] && atom[i].bo[j] != 9 )
320                                {
321                                        ib2 = atom[i].iat[j];
322                                        if( ib2 > ib1 )
323                                        {
324                                                an += 1.0F;
325                                                xx = atom[ib1].x - atom[ib2].x;
326                                                yy = atom[ib1].y - atom[ib2].y;
327                                                zz = atom[ib1].z - atom[ib2].z;
328                                                xyzdis += sqrt( xx*xx + yy*yy + zz*zz );
329                                        }
330                                }
331                        }
332                }
333                if (an > 0.0)
334                   flags.avleg = 1.10F/(xyzdis/an); /* no c-c bonds assume c-h bond dist */
335                else
336                   flags.avleg = 1.20F; /*  default when all else goes wrong */
337        }
338        else
339        {
340                flags.avleg = 1.20F; /*  default when all else goes wrong */
341        }
342        return;
343 }
344 /* -------------------------- */
345 void center(int istruct)
346 {
347        //  istruct = -2 do all
348        //  istruct = 1 do main structrure
349        //  istruct = number do that substructure
350
351        float centerx, centery, centerz;
352        int i;
353
354        centerx = 0.0F;
355        centery = 0.0F;
356        centerz = 0.0F;
357
358        if (istruct == -2)
359        {
360                for (i=1; i <= natom; i++)
361                {
362              centerx += atom[i].x/natom;
363          centery += atom[i].y/natom;
364                  centerz += atom[i].z/natom;
365                }
366                for (i=1; i <= natom; i++)
367                {
368                        atom[i].x -= centerx;
369                        atom[i].y -= centery;
370                        atom[i].z -= centerz;
371                }
372        } else if (istruct == -1)      // main structure only
373        {
374                for (i=1; i <= natom; i++)
375                {
376                        if(atom[i].substr[0] == 0)
377                        {
378                      centerx += atom[i].x/natom;
379                  centery += atom[i].y/natom;
380                          centerz += atom[i].z/natom;
381                        }
382                }
383                for (i=1; i <= natom; i++)
384                {
385                        if(atom[i].substr[0] == 0)
386                        {
387                                atom[i].x -= centerx;
388                                atom[i].y -= centery;
389                                atom[i].z -= centerz;
390                        }
391                }
392        } else                  // substructure
393        {
394                for (i=1; i <= natom; i++)
395                {
396                        if (atom[i].substr[0] & (1L << istruct) )
397                        {
398                      centerx += atom[i].x/natom;
399                  centery += atom[i].y/natom;
400                          centerz += atom[i].z/natom;
401                        }
402                }
403                for (i=1; i <= natom; i++)
404                {
405                        if (atom[i].substr[0] & (1L << istruct) )
406                        {
407                                atom[i].x -= centerx;
408                                atom[i].y -= centery;
409                                atom[i].z -= centerz;
410                        }
411                }
412        }
413 }
414 /* ------------------------ */
415
416 double ooplane(int ia,int ib,int ic,int id)
417 {
418      double angle,dot,cosine;
419      double xia,yia,zia;
420      double xib,yib,zib;
421      double xic,yic,zic;
422      double xid,yid,zid;
423      double xab,yab,zab,rab2,rab;
424      double xbc,ybc,zbc,rbc2,rbc;
425      double xbd,ybd,zbd,rbd2,rbd;
426      double sine;
427
428        /****the arguments of this function are the atom indices of atoms a-b-c-d
429         *     which determine dihedral angle dihdrl.
430         *     dihdrl returns the degree value of the dihedral('torsional')
431         *     angle defined by atoms i1,i2,i3, &i4.
432         */
433
434                 xia = atom[ia].x;
435                 yia = atom[ia].y;
436                 zia = atom[ia].z;
437                 xib = atom[ib].x;
438                 yib = atom[ib].y;
439                 zib = atom[ib].z;
440                 xic = atom[ic].x;
441                 yic = atom[ic].y;
442                 zic = atom[ic].z;
443                 xid = atom[id].x;
444                 yid = atom[id].y;
445                 zid = atom[id].z;
446                
447                 xab = xia - xib;
448                 yab = yia - yib;
449                 zab = zia - zib;
450                 rab2 = (xab*xab+yab*yab+zab*zab);
451                 rab = sqrt(rab2);
452                
453                 xbc = xic - xib;
454                 ybc = yic - yib;
455                 zbc = zic - zib;
456                 rbc2 = xbc*xbc+ybc*ybc+zbc*zbc;
457                 rbc = sqrt(rbc2);
458                
459                 if (rab2 != 0.0 && rbc2 != 0.0)
460                 {
461                    dot = xab*xbc + yab*ybc + zab*zbc;
462                    cosine = dot /(rab*rbc);
463                    if (cosine < -1.0)
464                       cosine = -1.0;
465                    if (cosine > 1.0)
466                       cosine = 1.0;
467                    sine = sqrt(1.00-cosine*cosine);
468
469                    xbd = atom[id].x - atom[ib].x;
470                    ybd = atom[id].y - atom[ib].y;
471                    zbd = atom[id].z - atom[ib].z;
472                    rbd2 = xbd*xbd+ybd*ybd+zbd*zbd;
473                    rbd = sqrt(rbd2);
474                    
475                    cosine = ( xbd*(yab*zbc-zab*ybc) + ybd*(zab*xbc-xab*zbc) + zbd*(xab*ybc-yab*xbc))/((rab*rbc*rbd)*sine);
476
477                    angle = radian * asin(cosine);
478              }
479        return (angle);
480 }
255   /* ------------------------------------------------------------- */      
256   void quick_type()
257   {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines