ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kcoord.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (11 years, 10 months ago) by tjod
File size: 3941 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2     #include "pcwin.h"
3     #include "pcmod.h"
4     #include "field.h"
5    
6     EXTERN struct t_vdw1 {
7     int nvdw;
8     float rad[MAXVDWCONST], eps[MAXVDWCONST];
9     int lpd[MAXVDWCONST], ihtyp[MAXVDWCONST], ihdon[MAXVDWCONST];
10     float alpha[MAXVDWCONST],n[MAXVDWCONST],a[MAXVDWCONST],g[MAXVDWCONST];
11     char da[MAXVDWCONST][2];
12     } vdw1;
13    
14     EXTERN struct t_cbonds {
15     int nbnd, icbonds[MAXCBND][3], lpd[MAXCBND], ia1[MAXCBND],ia2[MAXCBND];
16     double vrad[MAXCBND],veps[MAXCBND];
17     } cbonds;
18    
19     EXTERN int Missing_constants;
20     EXTERN FILE *errfile;
21    
22     void kcoord_bonds()
23     {
24     int i, ia, ib,ierr, ita,itb;
25     int jj, nusi;
26     long int mask5,mask6,mask7,mask8;
27     float rad,eps;
28    
29     mask5 = 1L << SATMET_MASK;
30     mask6 = 1L << GT18e_MASK;
31     mask7 = 1L << LOWSPIN_MASK;
32     mask8 = 1L << SQPLAN_MASK;
33    
34     for (i=0; i < cbonds.nbnd; i++)
35     {
36     ia = cbonds.icbonds[i][0];
37     ib = cbonds.icbonds[i][1];
38     ita = atom[ia].type;
39     itb = atom[ib].type;
40     nusi = 0;
41     rad = vdw1.rad[ita];
42     eps = vdw1.eps[ita];
43     ierr = FALSE;
44     if (strcmp(field.radiusrule,"ARITHMETIC") == 0)
45     {
46     cbonds.vrad[i] = rad + vdw1.rad[atom[ib].type];
47     } else if (strcmp(field.radiusrule,"GEOMETRIC") == 0)
48     {
49     cbonds.vrad[i] = 2.0*( sqrt(rad*vdw1.rad[atom[ib].type]));
50     } else
51     cbonds.vrad[i] = rad + vdw1.rad[atom[ib].type];
52    
53     if (strcmp(field.epsrule,"ARITHMETIC") == 0)
54     {
55     cbonds.veps[i] = 0.5*(eps + vdw1.eps[atom[ib].type]);
56     } else if (strcmp(field.epsrule,"GEOMETRIC") == 0)
57     {
58     cbonds.veps[i] = sqrt( eps*vdw1.eps[atom[ib].type] );
59     }else
60     cbonds.veps[i] = sqrt( eps*vdw1.eps[atom[ib].type] );
61    
62     if ( (atom[ia].flags & mask6) && !( atom[ia].flags & mask5 ))
63     nusi = 2;
64     else if( ( atom[ia].flags &mask5 ) && ( atom[ia].flags & mask6 ) )
65     {
66     if( ( atom[ia].flags & mask7 ) && !( atom[ia].flags & mask8 ) )
67     nusi = 1;
68     else if(( atom[ia].flags & mask8 ) && !( atom[ia].flags & mask7 ) )
69     nusi = 3;
70     }
71    
72     // metal to pi atom
73     if (atom[ib].type == 2 || atom[ib].type == 4 || atom[ib].type == 29 ||
74     atom[ib].type == 30 || atom[ib].type == 40 || atom[ib].type == 48 )
75     {
76     cbonds.ia1[i] = 0;
77     cbonds.ia2[i] = 0;
78     for (jj = 0 ; jj < MAXIAT; jj++)
79     {
80     if (atom[ib].iat[jj] != 0 && atom[ib].iat[jj] != ia && cbonds.ia1[i] == 0)
81     cbonds.ia1[i] = atom[ib].iat[jj];
82     else if (atom[ib].iat[jj] != 0 && atom[ib].iat[jj] != ia && cbonds.ia2[i] == 0)
83     cbonds.ia2[i] = atom[ib].iat[jj];
84     }
85     cbonds.vrad[i] = rad + 0.66;
86     if (nusi > 0)
87     cbonds.vrad[i] += 0.15;
88     cbonds.lpd[i] = vdw1.lpd[atom[ib].type];
89     if (nusi > 0)
90     cbonds.lpd[i] *= 0.85;
91     } else
92     { // metal to nitrogen and oxygen lone pairs
93     cbonds.vrad[i] = rad + vdw1.rad[atom[ib].type]; // use radius of lone pair
94     cbonds.vrad[i] *= 0.6; // cut radius in half
95    
96     cbonds.lpd[i] = vdw1.lpd[atom[ib].type];
97     if (atom[ib].type == 7)
98     cbonds.lpd[i] *= 3.0;
99     else if (atom[ib].type == 10)
100     cbonds.lpd[i] *= 4.0;
101    
102     if( ( atom[ia].flags & mask6 ) && !( atom[ia].flags & mask5 ) )
103     {
104     cbonds.vrad[i] += .15;
105     cbonds.lpd[i] *= .85;
106     }
107    
108     if( nusi == 1 || (nusi == 3 && (atom[atom[ib].iat[0]].type == 6 || atom[atom[ib].iat[0]].type == 15)) )
109     cbonds.vrad[i] -= .1;
110     }
111     }
112    
113     }