1 #include <chem/vibanal.h>
2 #include <madness/tensor/tensor_lapack.h>
3
4 using namespace madness;
5
6 /// returns the vibrational frequencies
7
8 /// @param[in] hessian the hessian matrix (not mass-weighted)
9 /// @param[out] normalmodes the normal modes
10 /// @param[in] project_tr whether to project out translation and rotation
11 /// @param[in] print_hessian whether to print the hessian matrix
12 /// @return the frequencies in atomic units
compute_frequencies(const Molecule & molecule,const Tensor<double> & hessian,Tensor<double> & normalmodes,const bool project_tr,const bool print_hessian)13 Tensor<double> compute_frequencies(const Molecule& molecule,
14 const Tensor<double>& hessian, Tensor<double>& normalmodes,
15 const bool project_tr, const bool print_hessian) {
16
17 // compute mass-weighing matrices
18 Tensor<double> M=molecule.massweights();
19 Tensor<double> Minv(3*molecule.natom(),3*molecule.natom());
20 for (int i=0; i<3*molecule.natom(); ++i) Minv(i,i)=1.0/M(i,i);
21
22 // mass-weight the hessian
23 Tensor<double> mwhessian=inner(M,inner(hessian,M));
24
25 // remove translation and rotation
26 if (project_tr) remove_external_dof(mwhessian,molecule);
27
28 if (print_hessian) {
29 if (project_tr) {
30 print("mass-weighted hessian with translation and rotation projected out");
31 } else {
32 print("mass-weighted unprojected hessian");
33 }
34 Tensor<double> mmhessian=inner(Minv,inner(mwhessian,Minv));
35 print(mwhessian);
36 print("mass-weighted unprojected hessian; mass-weighing undone");
37 print(mmhessian);
38 }
39
40 Tensor<double> freq;
41 syev(mwhessian,normalmodes,freq);
42 for (long i=0; i<freq.size(); ++i) {
43 if (freq(i)>0.0) freq(i)=sqrt(freq(i)); // real frequencies
44 else freq(i)=-sqrt(-freq(i)); // imaginary frequencies
45 }
46 return freq;
47 }
48
49
compute_reduced_mass(const Molecule & molecule,const Tensor<double> & normalmodes)50 Tensor<double> compute_reduced_mass(const Molecule& molecule, const Tensor<double>& normalmodes) {
51
52 Tensor<double> M=molecule.massweights();
53 Tensor<double> D=projector_external_dof(molecule);
54 Tensor<double> L=copy(normalmodes);
55 Tensor<double> DL=inner(D,L);
56 Tensor<double> MDL=inner(M,DL);
57 Tensor<double> mu(3*molecule.natom());
58
59 for (int i=0; i<3*molecule.natom(); ++i) {
60 double mu1=0.0;
61 for (int j=0; j<3*molecule.natom(); ++j) mu1+=MDL(j,i)*MDL(j,i);
62 if (mu1>1.e-14) mu(i)=1.0/(mu1*constants::atomic_mass_in_au);
63 }
64 return mu;
65 }
66
67
68 /// compute the projector to remove transl. and rot. degrees of freedom
69
70 /// taken from http://www.gaussian.com/g_whitepap/vib.htm
71 /// I don't really understand the concept behind the projectors, but it
72 /// seems to work, and it is not written down explicitly anywhere.
73 /// NOTE THE ERROR IN THE FORMULAS ON THE WEBPAGE !
projector_external_dof(const Molecule & mol)74 Tensor<double> projector_external_dof(const Molecule& mol) {
75
76 // compute the translation vectors
77 Tensor<double> transx(3*mol.natom());
78 Tensor<double> transy(3*mol.natom());
79 Tensor<double> transz(3*mol.natom());
80 for (int i=0; i<mol.natom(); ++i) {
81 transx[3*i ]=sqrt(mol.get_atom(i).get_mass_in_au());
82 transy[3*i+1]=sqrt(mol.get_atom(i).get_mass_in_au());
83 transz[3*i+2]=sqrt(mol.get_atom(i).get_mass_in_au());
84 }
85
86 // compute the rotation vectors
87
88 // move the molecule to its center of mass and compute
89 // the moment of inertia tensor
90 Tensor<double> com=mol.center_of_mass();
91 Molecule mol2=mol;
92 mol2.translate(-1.0*com);
93 Tensor<double> I=mol2.moment_of_inertia();
94 I.scale(constants::atomic_mass_in_au);
95
96 // diagonalize the moment of inertia
97 Tensor<double> v,e;
98 syev(I, v, e); // v being the "X" tensor on the web site
99 v=transpose(v);
100
101 // Tensor<double> B(e.size());
102 // for (long i=0; i<e.size(); ++i) B(i)=1.0/(2.0*e(i));
103 // print("rotational constants in cm-1");
104 // print(constants::au2invcm*B);
105
106 // rotation vectors
107 Tensor<double> rotx(3*mol.natom());
108 Tensor<double> roty(3*mol.natom());
109 Tensor<double> rotz(3*mol.natom());
110
111 for (int iatom=0; iatom<mol.natom(); ++iatom) {
112
113 // coordinates wrt the center of mass ("R" on the web site)
114 Tensor<double> coord(3);
115 coord(0l)=mol.get_atom(iatom).x-com(0l);
116 coord(1l)=mol.get_atom(iatom).y-com(1l);
117 coord(2l)=mol.get_atom(iatom).z-com(2l);
118
119 // note the wrong formula on the Gaussian website:
120 // multiply with sqrt(mass), do not divide!
121 coord.scale(sqrt(mol.get_atom(iatom).get_mass_in_au()));
122
123 // p is the dot product of R and X on the web site
124 Tensor<double> p=inner(coord,v);
125
126 // Eq. (5)
127 rotx(3*iatom + 0)=p(1)*v(0,2)-p(2)*v(0,1);
128 rotx(3*iatom + 1)=p(1)*v(1,2)-p(2)*v(1,1);
129 rotx(3*iatom + 2)=p(1)*v(2,2)-p(2)*v(2,1);
130
131 roty(3*iatom + 0)=p(2)*v(0,0)-p(0l)*v(0,2);
132 roty(3*iatom + 1)=p(2)*v(1,0)-p(0l)*v(1,2);
133 roty(3*iatom + 2)=p(2)*v(2,0)-p(0l)*v(2,2);
134
135 rotz(3*iatom + 0)=p(0l)*v(0,1)-p(1)*v(0,0);
136 rotz(3*iatom + 1)=p(0l)*v(1,1)-p(1)*v(1,0);
137 rotz(3*iatom + 2)=p(0l)*v(2,1)-p(1)*v(2,0);
138
139 }
140
141 // move the translational and rotational vectors to a common tensor
142 Tensor<double> ext_dof(6,3*mol.natom());
143 ext_dof(0l,_)=transx;
144 ext_dof(1l,_)=transy;
145 ext_dof(2l,_)=transz;
146 ext_dof(3l,_)=rotx;
147 ext_dof(4l,_)=roty;
148 ext_dof(5l,_)=rotz;
149
150 // normalize
151 for (int i=0; i<6; ++i) {
152 double norm=ext_dof(i,_).normf();
153 if (norm>1.e-14) ext_dof(i,_).scale(1.0/norm);
154 else ext_dof(i,_)=0.0;
155 }
156
157 // compute overlap to orthonormalize the projectors
158 Tensor<double> ovlp=inner(ext_dof,ext_dof,1,1);
159 syev(ovlp,v,e);
160 ext_dof=inner(v,ext_dof,0,0);
161
162 // normalize or remove the dof if necessary (e.g. linear molecules)
163 for (int i=0; i<6; ++i) {
164 if (e(i)<1.e-14) {
165 ext_dof(i,_).scale(0.0); // take out this degree of freedom
166 } else {
167 ext_dof(i,_).scale(1.0/sqrt(e(i))); // normalize
168 }
169 }
170
171 // construct projector on the complement of the rotations
172 Tensor<double> projector(3*mol.natom(),3*mol.natom());
173 for (int i=0; i<3*mol.natom(); ++i) projector(i,i)=1.0;
174
175 // compute the outer products of the projectors
176 // 1- \sum_i | t_i >< t_i |
177 projector-=inner(ext_dof,ext_dof,0,0);
178
179 // construct random tensor for orthogonalization
180 Tensor<double> D(3*mol.natom(),3*mol.natom());
181 D.fillrandom();
182 for (int i=0; i<6; ++i) D(i,_)=ext_dof(i,_);
183
184 // this works only for non-linear molecules -- projector seems simpler
185 // ovlp=inner(D,D,1,1);
186 // cholesky(ovlp); // destroys ovlp
187 // Tensor<double> L=transpose(ovlp);
188 // Tensor<double> Linv=inverse(L);
189 // D=inner(Linv,D,1,0);
190 // ovlp=inner(D,D,1,1);
191 //
192 // for (int i=0; i<6; ++i) D(i,_)=0.0;
193 // D=copy(D(Slice(6,8),_));
194
195 // return transpose(D);
196
197 return projector;
198
199 }
200
201 /// remove translational degrees of freedom from the hessian
remove_external_dof(Tensor<double> & hessian,const Molecule & mol)202 void remove_external_dof(Tensor<double>& hessian,
203 const Molecule& mol) {
204
205 // compute the translation of the center of mass
206 Tensor<double> projector_ext=projector_external_dof(mol);
207
208 // this is P^T * H * P
209 hessian=inner(projector_ext,inner(hessian,projector_ext),0,0);
210 }
211