ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/freemol/branches/sync4pymol12/src/mengine/src/vibrate.c

#	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