/**************************************************************************** * * DFT++: density functional package developed by * the research group of Prof. Tomas Arias, MIT. * * Principal author: Sohrab Ismail-Beigi * * Modifications for MPI version: Kenneth P Esler, * Sohrab Ismail-Beigi, and * Tairan Wang. * * Modifications for LSD version: Jason A Cline, * Gabor Csanyi, * Tairan Wang, and * Evan Reed * * Modifications for lattice/Pulay forces: Gabor Csanyi and * Sohrab Ismail-Beigi * * K-point symmetries: Nicolaj Moll * * Ionic dynamics and relaxation: Tairan Wang and * Evan Reed * * Copyright (C) 1996-2000 Sohrab Ismail-Beigi * * Permission granted to the academic community to * freely distribute and use this software provided that this notice is * maintained intact, and all publications resulting from its use * cite the following reference: * * Sohrab Ismail-Beigi and T. A. Arias, "New Algebraic Formulation of * Density Functional Calculation", Computer Physics Communications, * volume 128, issue 1-2, pages 1-45, June 2000. * ****************************************************************************/ /* * Code by Gabor Csanyi, May 2001 * * paw projectors * */ /* $Id: projector.cpp_defunct,v 1.1.2.2 2003/03/31 22:09:26 sahak Exp $ */ #include "header.h" /* * This code that calculates the projectors * for a given atom and a given angular momentum channel * (supports multiple projectors) * * it returns the result in a column bundle, where the number * of columns is the number of projectors for the channel * (usually 1, but not always) * * the column is in G-space, ready to do proj^C (= ) * */ void PAW_projector(int sp,int ion,int ll, int m, QuantumNumber *qnum, PW_Basis *basis, Lattice *lattice, Ioninfo *ioninfo, ColumnBundle &proj) { complex S,Ylm; real dot,posx,posy,posz,kx,ky,kz,kplusGx,kplusGy,kplusGz; real qx,qy,qz,q,invq,qhatx,qhaty,qhatz; real flq=0.0,r,invdq=0.0; real GGTxx,GGTyy,GGTzz,GGTxy,GGTxz,GGTyz; real Vol,invsqrtVol; real ym1,y0,y1,y2,a,b,c,x; int np,p,n,j,nmax=0,l; matrix3 G; /* pi et al. */ const real pi = M_PI; const real twopi = (real)2.0*pi; const real fourpi = (real)4.0*pi; const real sixth = (real)1.0/(real)6.0; /* Ylm normalization constants */ const real N00 = sqrt(1.0/(4.0*pi)); const real N10 = sqrt(3.0/(4.0*pi)); const real N11 = -sqrt(3.0/(8.0*pi)); const real N1m1 = sqrt(3.0/(8.0*pi)); const real N20 = sqrt(5.0/(16.0*pi)); const real N21 = -sqrt(15.0/(8.0*pi)); const real N2m1 = sqrt(15.0/(8.0*pi)); const real N22 = sqrt(15.0/(32.0*pi)); const real N2m2 = sqrt(15.0/(32.0*pi)); Speciesinfo *sp_info = &(ioninfo->species[sp]); /* proj will contain the final result */ np = proj.tot_ncols; proj.zero_out(); /* Get ionic position, k-vector, cell volume, G, and GGT matrix entries */ posx = sp_info->atpos[ion].v[0]; posy = sp_info->atpos[ion].v[1]; posz = sp_info->atpos[ion].v[2]; kx = qnum->kvec.v[0]; ky = qnum->kvec.v[1]; kz = qnum->kvec.v[2]; G = lattice->G; GGTxx = lattice->GGT.m[0][0]; GGTyy = lattice->GGT.m[1][1]; GGTzz = lattice->GGT.m[2][2]; GGTxy = lattice->GGT.m[0][1]; GGTxz = lattice->GGT.m[0][2]; GGTyz = lattice->GGT.m[1][2]; Vol = lattice->unit_cell_volume; invsqrtVol = (real)1.0/sqrt(Vol); /* Get l,m values for the potential */ l = sp_info->p_l[ll]; if (l < 0 || l > 2) die("l != {0,1,2} projectors not supported yet!\n"); /* Calculate the non-local pseudopotential matrix entries */ for (p = 0; p < np; p++) { invdq = (real)1.0/sp_info->p_dq_nl[ll][p]; nmax = sp_info->p_ngrid_nl[ll][p]-3; for (n = 0; n < basis->nbasis; n++) { /* q = k+G: here q is in reciprocal lattice coordinates */ qx = kplusGx = kx + basis->Gx[n]; qy = kplusGy = ky + basis->Gy[n]; qz = kplusGz = kz + basis->Gz[n]; /* S = exp(-i*dot(q,pos)) */ dot = -twopi*(posx*qx + posy*qy + posz*qz); S.x = cos(dot); S.y = sin(dot); q = sqrt( qx*qx*GGTxx + qy*qy*GGTyy + qz*qz*GGTzz + 2.0*(qx*qy*GGTxy + qx*qz*GGTxz + qy*qz*GGTyz) ); /* Interpolate flq using internal tables */ r = q*invdq; j = (int)r; if (j > nmax) { dft_log(DFT_SILENCE, "G=[%d %d %d] q=[%g %g %g] k=[%g %g %g]\n", basis->Gx[n],basis->Gy[n],basis->Gz[n], qx,qy,qz,kx,ky,kz); die("projector grid is too small!\n"); } /* Do cubic interpolation on flq with j-1, j, j+1, and j+2 * points from internal tables. If j==0, we clearly can't * real j-1. However f_l(-q) = (-1)^l*f_l(q) so we can * fill in the j-1 entry with the j+1 entry times (-1)^l */ if (j == 0) ym1 = (1-2*(l%2))*sp_info->p_flq[ll][p][j+1]; else ym1 = sp_info->p_flq[ll][p][j-1]; y0 = sp_info->p_flq[ll][p][j ]; y1 = sp_info->p_flq[ll][p][j+1]; y2 = sp_info->p_flq[ll][p][j+2]; a = sixth*(y2-ym1+3.0*(y0-y1)); b = 0.5*(y1+ym1-y0-y0); c = y1-sixth*(y2+3.0*y0+ym1+ym1); x = r-j; flq = ((a*x+b)*x+c)*x+y0; /* The line below is quadratic interpolation: */ /* flq = y0 + 0.5*x*(y1-ym1+x*(y1+ym1-y0-y0)); */ /* If we're at q=0, then we must handle the case separately * since Ylm and qhat are not defined for q=0. The cases are: * if l==0, Ylm is const. * if l>0, flq=0 anyways, so even though Ylm is undefined, * it doesn't matter what it is since the final result is zero. */ if (q < (real)1.0e-13) { if (l==0) { Ylm.x = N00; Ylm.y = 0.0; } else { Ylm.x = Ylm.y = 0.0; } goto skipcalcYlm; } /* Replace the above if block by the one below to reporduce * cgrad results. */ /* if (j == 0) */ /* { */ /* flq = y0 + (y1-y0)*x; */ /* if (l==0) { Ylm.x = N00; Ylm.y = 0.0; } */ /* else { Ylm.x = Ylm.y = 0.0; } */ /* goto skipcalcYlm; */ /* } */ /* Otherwise, we have to calculate Ylm: * qhat is a unit vector in direction of q: here q is in "real" * bohr^-1 units and q = (k+G)*G where (k+G) is a row-vector. */ qx = kplusGx*G.m[0][0] + kplusGy*G.m[1][0] + kplusGz*G.m[2][0]; qy = kplusGx*G.m[0][1] + kplusGy*G.m[1][1] + kplusGz*G.m[2][1]; qz = kplusGx*G.m[0][2] + kplusGy*G.m[1][2] + kplusGz*G.m[2][2]; invq = (real)1.0/q; qhatx = qx*invq; qhaty = qy*invq; qhatz = qz*invq; /* Calculate Ylm */ if (l==0) { Ylm.x = N00; Ylm.y = 0.0; } else if (l==1) { if (m==0) { Ylm.x = N10*qhatz; Ylm.y = 0.0; } else if (m==1) { Ylm.x = N11*qhatx; Ylm.y = N11*qhaty; } else if (m==-1) { Ylm.x = N1m1*qhatx; Ylm.y = -N1m1*qhaty; } } else if (l==2) { if (m==0) { Ylm.x = N20*(3.0*qhatz*qhatz-1.0); Ylm.y = 0.0; } else if (m==1) { Ylm.x = N21*(qhatz*qhatx); Ylm.y = N21*(qhatz*qhaty); } else if (m==-1) { Ylm.x = N2m1*( qhatz*qhatx); Ylm.y = N2m1*(-qhatz*qhaty); } else if (m==2) { Ylm.x = N22*(qhatx*qhatx-qhaty*qhaty); Ylm.y = N22*(2.0*qhatx*qhaty); } else if (m==-2) { Ylm.x = N2m2*(qhatx*qhatx-qhaty*qhaty); Ylm.y = N2m2*(-2.0*qhatx*qhaty); } } else die("projector: l ?!?!!\n"); skipcalcYlm: ; /* Calculate Vnl */ proj.col[p].c[n] = (S*Ylm); proj.col[p].c[n] *= fourpi*flq*invsqrtVol; } } /* We know at least these two things...set them to be safe. */ proj.qnum = qnum; } /* * This is the main driver to calculate PAW projections * it loops over atoms and channels, and * calls the routine above to get the projector. * */ void calc_projection(Everything &everything) { Elecinfo &einfo = everything.elecinfo; Elecvars &evars = everything.elecvars; Ioninfo &ioninfo = everything.ioninfo; Basis *basis = everything.basis_spec.basis; Lattice &lattice = everything.basis_spec.lattice; Symmetries &symm = everything.symm; column_bundle *C = evars.C; column_bundle *Y = evars.Y; dft_log("\n>>>>>>>>>>>>>>> Projections <<<<<<<<<<<<<<<<\n\n"); dft_log_flush(); dft_log("Diagonalizing Hsub ..."); dft_log_flush(); calc_UVCn(einfo, evars, symm, lattice); calc_Hsub(everything); for (int q=0; q < einfo.nstates; q++){ do_column_bundle_matrix_mult(evars.C[q], evars.Hsub_evecs[q], evars.Y[q], 0); evars.C[q] = evars.Y[q]; } dft_log("done\n\n"); dft_log_flush(); /* * set up projector dump file * * this should really be dumped and stamped with the other files... * */ FILE *projfile; projfile = dft_fopen("PROJ", "w"); // loop over the states: k points and spin for (int q = 0; q < einfo.nstates; q++) { // loop over the various species for(int sp = 0; sp < ioninfo.nspecies; sp++){ if(ioninfo.species[sp].projectorp == FALSE) continue; // loop over the ions in this species for(int ion = 0; ion < ioninfo.species[sp].natoms; ion++){ dft_log("q=%3d,k=[%3f %3f %3f],S=%2d,sp=%2d,ion=%3d ... ", q, einfo.qnums[q].kvec.v[0], einfo.qnums[q].kvec.v[1], einfo.qnums[q].kvec.v[2], einfo.qnums[q].spin, sp, ion); dft_log_flush(); //count up the number of projectors: loop over l, m, np (multiple projectors provision) int lmnp = 0; for (int l = 0; l < ioninfo.species[sp].p_nl; l++) for(int m = -l; m <= l; m++) lmnp += ioninfo.species[sp].p_np[l]; // allocate space for the projector column_bundle proj(lmnp, C[q].col_length); copy_innards_column_bundle(&(C[q]),&proj); /* dft_log("\n"); */ // fill in the projector int lmnp_i = 0; for (int l = 0; l < ioninfo.species[sp].p_nl; l++) for(int m = -l; m <= l; m++){ // create a temporary projector just for this l,m column_bundle lmproj(ioninfo.species[sp].p_np[l], C[q].col_length); copy_innards_column_bundle(&(C[q]),&lmproj); // fill in the project elements for this l,m PAW_projector(sp,ion,l,m,&(einfo.qnums[q]), &basis[q],&lattice,&ioninfo,lmproj); // copy the elements into the bundled projector for(int np = 0; np < ioninfo.species[sp].p_np[l]; np++) proj.col[lmnp_i++] = lmproj.col[np]; } // do the projection!! matrix pdagC(lmnp, einfo.nbands); pdagC = proj^C[q]; // dump the matrix to the file (append) pdagC.write(projfile); /* if(q == 0){ */ /* dft_log("\n"); */ /* complex c; */ /* for(int i = 0; i < lmnp; i++){ */ /* for(int j = 0; j < einfo.nbands; j++){ */ /* c = pdagC(i,j); */ /* dft_log("%8.4f+%8.4fi", c.x, c.y); */ /* } */ /* dft_log("\n"); */ /* } */ /* dft_log("\n"); */ /* } */ dft_log("done\n"); } } } dft_fclose(projfile); dft_log("\n>>>>>>>>>>>>>>> Projections Done !! <<<<<<<<<<<<<<<<\n\n"); dft_log_flush(); }