Commit 52ba2abd by Francois Gygi

fixed calculation with lower part of matrix.

added calculation using symmetrized matrix.
added calculation of eigenvectors.


git-svn-id: http://qboxcode.org/svn/qb/trunk@1543 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 53094b88
...@@ -18,7 +18,7 @@ ...@@ -18,7 +18,7 @@
// Forces are read from Qbox output // Forces are read from Qbox output
// The Qbox output should correspond to a sequence of calculations // The Qbox output should correspond to a sequence of calculations
// using symmetric finite displacements for all atoms in the x,y,z directions // using symmetric finite displacements for all atoms in the x,y,z directions
// Displacements have an amplitude of +/- h // Displacements have an amplitude of +/- h
// //
// use: dynmat force.dat h Nat1 mass1 [Nat2 mass2] ... // use: dynmat force.dat h Nat1 mass1 [Nat2 mass2] ...
// input_file: force.dat: forces from Qbox XML output file (collected with grep) // input_file: force.dat: forces from Qbox XML output file (collected with grep)
...@@ -40,8 +40,8 @@ ...@@ -40,8 +40,8 @@
#include <vector> #include <vector>
using namespace std; using namespace std;
extern "C" extern "C"
void dsyev_(const char *jobz, const char *uplo, const int *n, double *a, void dsyev_(const char *jobz, const char *uplo, const int *n, double *a,
const int *lda, double *w, double *wrk, int *lwrk, int *info); const int *lda, double *w, double *wrk, int *lwrk, int *info);
int main(int argc, char **argv) int main(int argc, char **argv)
...@@ -84,16 +84,16 @@ int main(int argc, char **argv) ...@@ -84,16 +84,16 @@ int main(int argc, char **argv)
{ {
infile >> dum >> fx >> fy >> fz; infile >> dum >> fx >> fy >> fz;
while (infile.get() != '\n'); while (infile.get() != '\n');
f[ishift][ia1][idir1][ia2][0] = fx; f[ishift][ia1][idir1][ia2][0] = fx;
f[ishift][ia1][idir1][ia2][1] = fy; f[ishift][ia1][idir1][ia2][1] = fy;
f[ishift][ia1][idir1][ia2][2] = fz; f[ishift][ia1][idir1][ia2][2] = fz;
} }
} }
} }
} }
// compute matrix elements: centered finite difference // compute matrix elements: centered finite difference
// a[3*ia1+idir1][3*ia2+idir2] = // a[3*ia1+idir1][3*ia2+idir2] =
// ( f[1][ia1][idir1][ia2][idir2] - f[0][ia1][idir1][ia2][idir2] ) / ( 2 h ) // ( f[1][ia1][idir1][ia2][idir2] - f[0][ia1][idir1][ia2][idir2] ) / ( 2 h )
const int n = 3*nat; const int n = 3*nat;
...@@ -109,7 +109,7 @@ int main(int argc, char **argv) ...@@ -109,7 +109,7 @@ int main(int argc, char **argv)
{ {
int i = 3*ia1+idir1; int i = 3*ia1+idir1;
int j = 3*ia2+idir2; int j = 3*ia2+idir2;
double aij = ( f[1][ia1][idir1][ia2][idir2] - double aij = ( f[1][ia1][idir1][ia2][idir2] -
f[0][ia1][idir1][ia2][idir2] ) / ( 2.0 * h ); f[0][ia1][idir1][ia2][idir2] ) / ( 2.0 * h );
aij = aij / sqrt(mass[ia1]*mass[ia2]); aij = aij / sqrt(mass[ia1]*mass[ia2]);
a[j*n+i] = aij; a[j*n+i] = aij;
...@@ -148,6 +148,9 @@ int main(int argc, char **argv) ...@@ -148,6 +148,9 @@ int main(int argc, char **argv)
a[i*n+j] = aij; a[i*n+j] = aij;
} }
} }
// compute eigenvectors
jobz = 'v';
dsyev_(&jobz,&uplo,&n,&a[0],&n,&w_sym[0],&wrk[0],&lwrk,&info); dsyev_(&jobz,&uplo,&n,&a[0],&n,&w_sym[0],&wrk[0],&lwrk,&info);
assert(info==0); assert(info==0);
...@@ -155,7 +158,7 @@ int main(int argc, char **argv) ...@@ -155,7 +158,7 @@ int main(int argc, char **argv)
for ( int i = 0; i < n; i++ ) for ( int i = 0; i < n; i++ )
{ {
if ( w_upper[i] < 0.0 ) if ( w_upper[i] < 0.0 )
cout << setw(8) << (int) (Ha2cm1 * sqrt(-w_upper[i])) << " I"; cout << setw(8) << (int) (Ha2cm1 * sqrt(-w_upper[i])) << " I";
else else
cout << setw(8) << (int) (Ha2cm1 * sqrt(w_upper[i])); cout << setw(8) << (int) (Ha2cm1 * sqrt(w_upper[i]));
...@@ -170,4 +173,24 @@ int main(int argc, char **argv) ...@@ -170,4 +173,24 @@ int main(int argc, char **argv)
cout << setw(8) << (int) (Ha2cm1 * sqrt(w_sym[i])); cout << setw(8) << (int) (Ha2cm1 * sqrt(w_sym[i]));
cout << endl; cout << endl;
} }
ofstream vecfile("dynmat_eigvec.dat");
for ( int j = 0; j < n; j++ )
{
vecfile << "# mode " << j+1 << " frequency = ";
if ( w_sym[j] < 0.0 )
vecfile << setw(8) << (int) (Ha2cm1 * sqrt(-w_sym[j])) << " I";
else
vecfile << setw(8) << (int) (Ha2cm1 * sqrt(w_sym[j]));
vecfile << " cm-1" << endl;
vecfile << setprecision(8);
vecfile.setf(ios::fixed, ios::floatfield);
vecfile.setf(ios::right, ios::adjustfield);
for ( int i = 0; i < n; i+=3 )
vecfile << i/3+1 << " "
<< setw(12) << a[j*n+i] << " "
<< setw(12) << a[j*n+i+1] << " "
<< setw(12) << a[j*n+i+2] << endl;
}
} }
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment