ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/vibrate.c
(Generate patch)
# Line 2 | Line 2
2  
3   #include "pcwin.h"
4   #include "pcmod.h"
5 +
6 + #include "asnsym.h"
7   #include "utility.h"
8   #include "hess.h"
9   #include "field.h"
10   #include "pot.h"
11   #include "energies.h"
12   #include "bonds_ff.h"
13 + #include "vibrate.h"
14  
15   char *symname[] = {"A","A'","A''","Ag","Au",
16                  "A1","A2","A1'","A2'","A1''","A2''",
# Line 18 | Line 21
21                  "Eg","Eu","E1'","E1''","E2'",
22                  "E2''","E1g","E1u","E2g","E2u","T","Tg","Tu" };
23  
24 < EXTERN struct {
22 <    int ishape,numsym,chiral;
23 <    double xi,yi,zi;
24 <    }  symmetry;
25 < EXTERN struct t_dipolemom {
26 <        double total, xdipole, ydipole, zdipole;
27 <       }  dipolemom;
28 <      
29 < struct t_quadmom {
30 <        double xx,xy,xz,yx,yy,yz,zx,zy,zz;
31 <       } quadmom;
32 <      
33 < EXTERN struct t_vibdata {
34 <        char ptgrp[4];
35 <        float mom_ix,mom_iy,mom_iz;
36 <        float etot,htot,stot,gtot,cptot;
37 <       } vibdata;
38 <      
24 >              
25   EXTERN double hesscut;
26  
27 < FILE * fopen_path ( char * , char * , char * ) ;
28 < double distance(int,int);
29 < void eigenxyz(double *, double **);
30 < void hessian(int, double *,int *, int *, int *,double *);
31 < void tred2(double **,int , double *, double *);
32 < void tqli(double *,double *,int ,double **);
33 < void message_alert(char *, char *);
48 < void ptgrp(int,char *,double **);
49 < void assign_vibsym(char *,char *,int *,double **);
50 < int findh(double **,double **,int);
51 < void vibcharge_dipole(double **);
52 < void vibdipole_dipole(double **);
53 < int findi(double **, double **);
54 < void vibrate(void);
27 > // FILE * fopen_path ( char * , char * , char * ) ;
28 >
29 > static void tred2(double **,int , double *, double *);
30 > static void tqli(double *,double *,int ,double **);
31 > static void assign_vibsym(char *,char *,int *,double **);
32 > //void vibcharge_dipole(double **);
33 > //void vibdipole_dipole(double **);
34  
35   void vibrate()
36   {
# Line 219 | Line 198
198            assign_vibsym(ngrp,vsym,&nsym,acoord);
199            vibsym[i] = nsym;
200            
222    //      dipolemom.total = (dqx*dqx + dqy*dqy + dqz*dqz)*97.217;
223    //      dipole[i] = dipolemom.total;
201        }
202        
203   /*      icycle = 3*natom/5;
# Line 282 | Line 259
259        }
260        fclose(vibfile); */
261  
262 < /*      fprintf(pcmoutfile,"\n==========================================================\n");
263 <      fprintf(pcmoutfile,"Normal Mode Vibrational Analysis\n\n");
262 > /*      fprintf(pcmlogfile,"\n==========================================================\n");
263 >      fprintf(pcmlogfile,"Normal Mode Vibrational Analysis\n\n");
264        if (symmetry.ishape == 1)
265 <         fprintf(pcmoutfile,"Molecular Shape: LINEAR\n");
265 >         fprintf(pcmlogfile,"Molecular Shape: LINEAR\n");
266        else if (symmetry.ishape == 2)
267 <         fprintf(pcmoutfile,"Molecular Shape: ASYMMETRIC ROTOR\n");
267 >         fprintf(pcmlogfile,"Molecular Shape: ASYMMETRIC ROTOR\n");
268        else if (symmetry.ishape == 3)
269 <         fprintf(pcmoutfile,"Molecular Shape: PROLATE SYMMETRIC ROTOR\n");
269 >         fprintf(pcmlogfile,"Molecular Shape: PROLATE SYMMETRIC ROTOR\n");
270        else if (symmetry.ishape == 4)
271 <         fprintf(pcmoutfile,"Molecular Shape: OBLATE SYMMETRIC ROTOR\n");
271 >         fprintf(pcmlogfile,"Molecular Shape: OBLATE SYMMETRIC ROTOR\n");
272        else if (symmetry.ishape == 5)
273 <         fprintf(pcmoutfile,"Molecular Shape: SPHERICAL ROTOR\n");
273 >         fprintf(pcmlogfile,"Molecular Shape: SPHERICAL ROTOR\n");
274        else if (symmetry.ishape == 6)
275 <         fprintf(pcmoutfile,"Molecular Shape: PLANAR\n");
275 >         fprintf(pcmlogfile,"Molecular Shape: PLANAR\n");
276        
277 <      fprintf(pcmoutfile,"Symmetry Point Group: %s\n",ngrp);
278 <      fprintf(pcmoutfile,"Symmetry Number: %d\n",symmetry.numsym);
277 >      fprintf(pcmlogfile,"Symmetry Point Group: %s\n",ngrp);
278 >      fprintf(pcmlogfile,"Symmetry Number: %d\n",symmetry.numsym);
279        if (symmetry.chiral)
280 <          fprintf(pcmoutfile,"Chirality: Yes\n");
280 >          fprintf(pcmlogfile,"Chirality: Yes\n");
281        else
282 <          fprintf(pcmoutfile,"Chirality: No\n"); */
282 >          fprintf(pcmlogfile,"Chirality: No\n"); */
283  
284        if (symmetry.ishape == 1)
285           mm = maxvib - 5;
286        else
287           mm = maxvib - 6;
288  
289 < /*      fprintf(pcmoutfile,"==========================================================\n");
290 <      fprintf(pcmoutfile,"\nFundmental Normal Vibrational Frequencies\n");
289 > /*      fprintf(pcmlogfile,"==========================================================\n");
290 >      fprintf(pcmlogfile,"\nFundmental Normal Vibrational Frequencies\n");
291        for (i=0; i < mm; i++)
292 <         fprintf(pcmoutfile,"   %-3d            %8.3f      %s    %8.3f\n",i,eigen[i],symname[vibsym[i]],dipole[i]);
292 >         fprintf(pcmlogfile,"   %-3d            %8.3f      %s    %8.3f\n",i,eigen[i],symname[vibsym[i]],dipole[i]);
293          
294 <      fprintf(pcmoutfile,"\nMoments of Inertia\n");
295 <      fprintf(pcmoutfile," IX = %8.3f   IY = %8.3f   IZ = %8.3f\n",symmetry.xi,symmetry.yi,symmetry.zi); */
294 >      fprintf(pcmlogfile,"\nMoments of Inertia\n");
295 >      fprintf(pcmlogfile," IX = %8.3f   IY = %8.3f   IZ = %8.3f\n",symmetry.xi,symmetry.yi,symmetry.zi); */
296              
297        strcpy(vibdata.ptgrp,ngrp);
298        vibdata.mom_ix = symmetry.xi;
# Line 414 | Line 391
391        vibdata.gtot = gtot;
392        vibdata.cptot = cptot;
393   // print out
394 < /*      fprintf(pcmoutfile,"\nTemperatur: %8.3f K \n",temp);
395 <      fprintf(pcmoutfile,"\nZero Point Energy: %8.3f kcal/mol \n",ezpe);
396 <      fprintf(pcmoutfile,"\n                    Energy       Enthalpy       Entropy      Free Energy      Heat Capacity\n");
397 <      fprintf(pcmoutfile,"\n                    kcal/mol     kcal/mol       eu           kcal/mol          cal/mol/deg\n");
398 <      fprintf(pcmoutfile,"Translational: %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",etrans,htrans,strans,gtrans,cptrans);
399 <      fprintf(pcmoutfile,"Rotational:    %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",erot,hrot,srot,grot,cprot);
400 <      fprintf(pcmoutfile,"Vibrational:   %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",evib,hvib,svib,gvib,cpvib);
401 <      fprintf(pcmoutfile,"Mixing:        %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",emix,hmix,smix,gmix,cpmix);
402 <      fprintf(pcmoutfile,"Total:         %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",etot,htot,stot,gtot,cptot); */
394 > /*      fprintf(pcmlogfile,"\nTemperatur: %8.3f K \n",temp);
395 >      fprintf(pcmlogfile,"\nZero Point Energy: %8.3f kcal/mol \n",ezpe);
396 >      fprintf(pcmlogfile,"\n                    Energy       Enthalpy       Entropy      Free Energy      Heat Capacity\n");
397 >      fprintf(pcmlogfile,"\n                    kcal/mol     kcal/mol       eu           kcal/mol          cal/mol/deg\n");
398 >      fprintf(pcmlogfile,"Translational: %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",etrans,htrans,strans,gtrans,cptrans);
399 >      fprintf(pcmlogfile,"Rotational:    %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",erot,hrot,srot,grot,cprot);
400 >      fprintf(pcmlogfile,"Vibrational:   %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",evib,hvib,svib,gvib,cpvib);
401 >      fprintf(pcmlogfile,"Mixing:        %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",emix,hmix,smix,gmix,cpmix);
402 >      fprintf(pcmlogfile,"Total:         %10.3f     %10.3f     %10.3f    %10.3f       %10.3f\n",etot,htot,stot,gtot,cptot); */
403        
404   //      
405      free_dvector(eigen, 0, 3*natom);
# Line 1137 | Line 1114
1114      }
1115      free_dmatrix(b ,0, natom+1, 0,3);
1116   }    
1117 + // ======================
1118 + void tred2(double **a,int n, double *d, double *e)
1119 + {
1120 +
1121 +    int i,j,k,l, nn;
1122 +    double scale, hh, h, g, f;
1123 +
1124 +    nn = n-1;
1125 +    for (i= nn; i >= 1; i--)
1126 +    {
1127 +        l = i - 1;
1128 +        h = scale = 0.0;
1129 +        if (l > 0)
1130 +        {
1131 +            for (k=0; k <= l; k++)
1132 +                scale += fabs(a[i][k]);
1133 +            if (scale == 0.0)
1134 +                e[i] = a[i][l];
1135 +            else
1136 +            {
1137 +                for (k=0; k <= l; k++)
1138 +                {
1139 +                    a[i][k] /= scale;
1140 +                    h += a[i][k]*a[i][k];
1141 +                }
1142 +                f = a[i][l];
1143 +                g = f>0 ? -sqrt(h) : sqrt(h);
1144 +                e[i] = scale*g;
1145 +                h -= f*g;
1146 +                a[i][l] = f-g;
1147 +                f = 0.0;
1148 +                for (j=0; j <= l; j++)
1149 +                {
1150 +                    a[j][i] = a[i][j]/h;
1151 +                    g = 0.0;
1152 +                    for (k=0; k <= j; k++)
1153 +                        g += a[j][k]*a[i][k];
1154 +                    for (k = j+1; k <= l; k++)
1155 +                        g += a[k][j]*a[i][k];
1156 +                    e[j] = g/h;
1157 +                    f += e[j]*a[i][j];
1158 +                }
1159 +                hh = f/(h+h);
1160 +                for (j=0; j <= l; j++)
1161 +                {
1162 +                    f = a[i][j];
1163 +                    e[j] = g = e[j] - hh*f;
1164 +                    for (k=0; k <= j; k++)
1165 +                      a[j][k] -= (f*e[k]+g*a[i][k]);
1166 +                }
1167 +            }
1168 +        }else
1169 +           e[i] = a[i][l];
1170 +        d[i] = h;
1171 +    }
1172 +    d[0] = 0.0;
1173 +    e[0] = 0.0;
1174 +    for (i=0; i < n; i++)
1175 +    {
1176 +        l = i - 1;
1177 +        if (i > 0)
1178 +        {
1179 +            for (j=0; j <= l; j++)
1180 +            {
1181 +                g = 0.0;
1182 +                for (k=0; k <= l; k++)
1183 +                    g += a[i][k]*a[k][j];
1184 +                for (k=0; k <= l; k++)
1185 +                    a[k][j] -= g*a[k][i];
1186 +            }
1187 +        }
1188 +        d[i] = a[i][i];
1189 +        a[i][i] = 1.0;
1190 +        for (j=0; j <= l; j++)
1191 +        {
1192 +            a[j][i] = a[i][j] = 0.0;
1193 +        }
1194 +    }
1195 + }
1196 + /* =========================================== */
1197 + void tqli(double *d, double *e, int n, double **z)
1198 + {
1199 +    int m,l,iter,i,k,j;
1200 +    double s,r,p,g,f,dd,c,b;
1201 +
1202 +    for (i=1; i < n; i++)
1203 +       e[i-1] = e[i];
1204 +    e[n-1] = 0.0;
1205 +
1206 +    for (l=0; l < n; l++)
1207 +    {
1208 +        iter = 0;
1209 +        do
1210 +        {
1211 +            for (m=l; m < n-1; m++)
1212 +            {
1213 +                dd = fabs(d[m]) + fabs(d[m+1]);
1214 +                if (fabs(e[m])+dd == dd) break;
1215 +            }
1216 +            if (m!= l)
1217 +            {
1218 +                if (iter++ == 30)
1219 +                {
1220 +                    fprintf(pcmlogfile,"Error in tqli\n");
1221 +                    message_alert("Error in Tqli","Error");
1222 +                    return;
1223 +                }
1224 +                g = (d[l+1]-d[l])/(2.0*e[l]);
1225 +                r = sqrt((g*g)+1.0);
1226 +                g = d[m]-d[l]+e[l]/(g+SIGN(r,g));
1227 +                s = c = 1.0;
1228 +                p = 0.0;
1229 +                for (i=m-1; i >= l; i--)
1230 +                {
1231 +                    f = s*e[i];
1232 +                    b = c*e[i];
1233 +                    if (fabs(f) >= fabs(g))
1234 +                    {
1235 +                        c = g/f;
1236 +                        r = sqrt((c*c)+1.0);
1237 +                        e[i+1] = f*r;
1238 +                        c *= (s=1.0/r);
1239 +                    } else
1240 +                    {
1241 +                        s = f/g;
1242 +                        r = sqrt((s*s)+1.0);
1243 +                        e[i+1] = g*r;
1244 +                        s *= (c=1.0/r);
1245 +                    }
1246 +                    g = d[i+1]-p;
1247 +                    r = (d[i]-g)*s+2.0*c*b;
1248 +                    p = s*r;
1249 +                    d[i+1] = g+p;
1250 +                    g = c*r-b;
1251 +                    for (k=0; k < n; k++)
1252 +                    {
1253 +                        f = z[k][i+1];
1254 +                        z[k][i+1] = s*z[k][i]+c*f;
1255 +                        z[k][i] = c*z[k][i]-s*f;
1256 +                    }
1257 +                }
1258 +                d[l] = d[l]-p;
1259 +                e[l] = g;
1260 +                e[m] = 0.0;
1261 +            }
1262 +        } while (m!= l);
1263 +    }
1264 + // sort d and z
1265 +    for (i=0; i < n; i++)
1266 +    {
1267 +        p  = d[k=i];
1268 +        for (j=i+1; j < n; j++)
1269 +            if (d[j] >= p) p = d[k=j];
1270 +        if (k != i)
1271 +        {
1272 +            d[k] = d[i];
1273 +            d[i] = p;
1274 +            for (j=0; j < n; j++)
1275 +            {
1276 +                p = z[j][i];
1277 +                z[j][i] = z[j][k];
1278 +                z[j][k] = p;
1279 +            }
1280 +        }
1281 +    }  
1282 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines