API tools faq
paste
ETikhonov

Calculation program

ETikhonov
Jun 17th, 2014
180
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 21.95 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <fstream>
  3. #include <list>
  4. #include <string>
  5. #include <sstream>
  6. #include <cmath>
  7. using std::fstream;
  8. using std::ofstream;
  9. using std::cout;
  10. using std::endl;
  11. using std::list;
  12. using std::string;
  13. using std::stringstream;
  14. double Integrate3(double *func, int Nx, int Ny, int Nz);
  15. double Integrate6(double *func, int Nx, int Ny, int Nz);
  16. void Renorm(double *func, int Nx, int Ny, int Nz, double norm);
  17. void SquareMO(double *source,double *dest, int Nx, int Ny, int Nz);
  18. void Substract(double *term1, double *term2, double *source,
  19.                int Nx, int Ny, int Nz);
  20. void Mult(double *term1, double *term2, double *source,
  21.                int Nx, int Ny, int Nz);
  22. void InitDeltaDens(double *source,double *dest, int Nx, int Ny, int Nz);
  23. void InitLeftDens(double *source,double *dest, double *diff,
  24.                   int Nx, int Ny, int Nz);
  25. void InitRightDens(double *source,double *dest, double *diff,
  26.                   int Nx, int Ny, int Nz);
  27. void InitRS(double *source,double *dest,
  28.                   int Nx, int Ny, int Nz);
  29. void InitEx(double *rho, double *rs, double *dest,
  30.                   int Nx, int Ny, int Nz);
  31. void InitLeftEx_m(double *rho, double *rs, double *dest,
  32.                   int Nx, int Ny, int Nz);
  33. void InitRightEx_m(double *rho, double *rs, double *dest,
  34.                   int Nx, int Ny, int Nz);
  35. void InitLeftEx_r(double *rho, double *rs, double *dest,
  36.                   int Nx, int Ny, int Nz);
  37. void InitRightEx_r(double *rho, double *rs, double *dest,
  38.                   int Nx, int Ny, int Nz);
  39. void InitEc(double *rho, double *rs, double *dest,
  40.                   int Nx, int Ny, int Nz);
  41. void InitLeftEc_m(double *rho, double *rs, double *dest,
  42.                   int Nx, int Ny, int Nz);
  43. void InitRightEc_m(double *rho, double *rs, double *dest,
  44.                   int Nx, int Ny, int Nz);
  45. void InitLeftEc_r(double *rho, double *rs, double *dest,
  46.                   int Nx, int Ny, int Nz);
  47. void InitRightEc_r(double *rho, double *rs, double *dest,
  48.                   int Nx, int Ny, int Nz);
  49. void FormIntM(double *rhoHOMO, double *rhoLUMO, double *rho,
  50.               double *Ex, double *Exleft, double *Exright,
  51.               double *Ec, double *Ecleft, double *Ecright,double *dest,
  52.               int Nx, int Ny, int Nz);
  53. void FormIntR(double *rhoHOMO, double *rhoLUMO, double *drho,
  54.               double *Ex, double *Exleft, double *Exright,
  55.               double *Ec, double *Ecleft, double *Ecright,double *dest,
  56.               int Nx, int Ny, int Nz);
  57. void WriteGCube(string filename, double *data,int Nx, int Ny, int Nz);
  58. void ConstructV(double *charge, double *dest, int Nx, int Ny, int Nz);
  59. double dx,dy,dz;
  60. double h=0.001; /* increment for differentiating f(z) */
  61. const double PI=3.141592653589793;
  62. /* PWLDA parameters below */
  63. double f_z=0.0;
  64. double f_lz=(pow(1-h,4/3.)+pow(1+h,4/3.)-2)/(pow(2,4/3.)-2);
  65. double f_rz=(pow(1+h,4/3.)+pow(1-h,4/3.)-2)/(pow(2,4/3.)-2);
  66. double f2z=1.709921;
  67. double AEX=-(3./(4*PI))*pow((9*PI/4.),1/3.);
  68. /* for \epsilon_c(r_s,0) */
  69. double A_0=0.031091;
  70. double alpha1_0=0.21370;
  71. double beta1_0=7.5957;
  72. double beta2_0=3.5876;
  73. double beta3_0=1.6382;
  74. double beta4_0=0.49294;
  75. double p=1.00; /* Common for all */
  76. /* for \epsilon_c(r_s,1) */
  77. double A_1=0.015545;
  78. double alpha1_1=0.20548;
  79. double beta1_1=14.1189;
  80. double beta2_1=6.1977;
  81. double beta3_1=3.3662;
  82. double beta4_1=0.62517;
  83. /* for \alpha_c(rs) */
  84. double A_a=0.016887;
  85. double alpha1_a=0.11125;
  86. double beta1_a=10.357;
  87. double beta2_a=3.6231;
  88. double beta3_a=0.88026;
  89. double beta4_a=0.49671;
  90. /* For writing CUB file these need to be global */
  91. int Nat=0;
  92. int Nx,Ny,Nz;
  93. double xmin,ymin,zmin;
  94. /* START OF main() */
  95. int main (int ac, char* av[]){
  96. if (ac != 4)
  97.     {
  98.         cout<<"Usage: "<<av[0]<< " <density>.cub <homo>.cub <lumo>.cub"<<endl;
  99.         return 1;
  100.     }
  101. cout.precision(16);
  102. fstream fdensity;
  103. fstream fhomo;
  104. fstream flumo;
  105. fdensity.open(av[1]);
  106. if (!fdensity.is_open())
  107. {
  108. cout<<"Cannot open file "<<av[1]<<"."<<endl;
  109. return 1;
  110. }
  111. cout << "Reading density from file "<<av[1]<<"... ";
  112. string buffer;
  113. int nline=0;
  114. while(getline(fdensity,buffer))
  115. {
  116. double trash;
  117. if (nline==2)
  118. {
  119. stringstream s(buffer);
  120. s>>Nat>>xmin>>ymin>>zmin;
  121. if (Nat<0) Nat=-Nat;
  122. }
  123. if (nline==3){
  124. stringstream s(buffer);
  125. s>>Nx>>dx>>trash>>trash;
  126. }
  127. if (nline==4){
  128. stringstream s(buffer);
  129. s>>Ny>>trash>>dy>>trash;
  130. }
  131. if (nline==5){
  132. stringstream s(buffer);
  133. s>>Nz>>trash>>trash>>dz;
  134. }
  135. if (nline>=5+Nat)
  136. {
  137. break;
  138. }
  139. nline++;
  140. }
  141. double *density = new double[Nx*Ny*Nz];
  142. int i,j,k;
  143. for(i=0;i<Nx;i++){
  144.   for(j=0;j<Ny;j++){
  145.     for(k=0;k<Nz;k++){
  146.       fdensity>>density[i*Ny*Nz+j*Nz+k];
  147. }}}
  148. fdensity.close();
  149. cout<<"Done!"<<endl;
  150. cout<<"Number of atoms = "<<Nat<<", grid points x = "<<Nx
  151.     <<", grid points y = "<<Ny<<", grid points z = "<<Nz<<endl;
  152. fhomo.open(av[2]);
  153. if (!fhomo.is_open())
  154. {
  155. cout<<"Cannot open file "<<av[2]<<"."<<endl;
  156. return 1;
  157. }
  158. cout << "Reading HOMO from file "<<av[2]<<"... ";
  159. nline=0;
  160. while(getline(fhomo,buffer)){
  161. if (nline==Nat+5) break;
  162. nline++;
  163. }
  164. double *homo = new double[Nx*Ny*Nz];
  165. double *homo2 = new double[Nx*Ny*Nz];
  166. for(i=0;i<Nx;i++){
  167.   for(j=0;j<Ny;j++){
  168.     for(k=0;k<Nz;k++){
  169.       fhomo>>homo[i*Ny*Nz+j*Nz+k];
  170. }}}
  171. fhomo.close();
  172. SquareMO(homo,homo2,Nx,Ny,Nz);
  173. cout<<"Done!"<<endl;
  174. flumo.open(av[3]);
  175. if (!flumo.is_open())
  176. {
  177. cout<<"Cannot open file "<<av[3]<<"."<<endl;
  178. return 1;
  179. }
  180. cout << "Reading LUMO from file "<<av[3]<<"... ";
  181. nline=0;
  182. while(getline(flumo,buffer)){
  183. if (nline==Nat+5) break;
  184. nline++;
  185. }
  186. double *lumo = new double [Nx*Ny*Nz];
  187. double *lumo2 = new double [Nx*Ny*Nz];
  188. for(i=0;i<Nx;i++){
  189.   for(j=0;j<Ny;j++){
  190.     for(k=0;k<Nz;k++){
  191.       flumo>>lumo[i*Ny*Nz+j*Nz+k];
  192. }}}
  193. flumo.close();
  194. //END of reading data
  195. SquareMO(lumo,lumo2,Nx,Ny,Nz);
  196. cout<<"Done!"<<endl;
  197. double HOMO_NORM=Integrate3(homo2,Nx,Ny,Nz);
  198. double LUMO_NORM=Integrate3(lumo2,Nx,Ny,Nz);
  199. cout<<"Integrated number of electrons = "
  200.     <<Integrate3(density,Nx,Ny,Nz)<<endl;
  201. cout<<"Integrated HOMO norm = "<<HOMO_NORM<<endl;
  202. cout<<"Integrated LUMO norm = "<<LUMO_NORM<<endl;
  203. cout<<"Renorming HOMO... ";
  204. Renorm(homo,Nx,Ny,Nz,sqrt(HOMO_NORM));
  205. SquareMO(homo,homo2,Nx,Ny,Nz);
  206. cout<<"new norm = "<<Integrate3(homo2,Nx,Ny,Nz)<<endl;
  207. cout<<"Renorming LUMO... ";
  208. Renorm(lumo,Nx,Ny,Nz,sqrt(LUMO_NORM));
  209. SquareMO(lumo,lumo2,Nx,Ny,Nz);
  210. cout<<"new norm = "<<Integrate3(lumo2,Nx,Ny,Nz)<<endl;
  211. double *delta_dens = new double [Nx*Ny*Nz];
  212. InitDeltaDens(density,delta_dens,Nx,Ny,Nz);
  213. double *left_dens = new double [Nx*Ny*Nz];
  214. InitLeftDens(density,left_dens,delta_dens,Nx,Ny,Nz);
  215. double *right_dens = new double [Nx*Ny*Nz];
  216. InitRightDens(density,right_dens,delta_dens,Nx,Ny,Nz);
  217. double *rs = new double [Nx*Ny*Nz];
  218. InitRS(density,rs,Nx,Ny,Nz);
  219. double *rsleft = new double [Nx*Ny*Nz];
  220. InitRS(left_dens,rsleft,Nx,Ny,Nz);
  221. double *rsright = new double [Nx*Ny*Nz];
  222. InitRS(right_dens,rsright,Nx,Ny,Nz);
  223. double *Ex = new double [Nx*Ny*Nz];
  224. InitEx(density,rs,Ex,Nx,Ny,Nz);
  225. cout<<"Exchange energy: "<<Integrate3(Ex,Nx,Ny,Nz)<<" Hartree"<<endl;
  226. double *left_Ex_m = new double [Nx*Ny*Nz];
  227. InitLeftEx_m(density,rs,left_Ex_m,Nx,Ny,Nz);
  228. double *right_Ex_m = new double [Nx*Ny*Nz];
  229. InitRightEx_m(density,rs,right_Ex_m,Nx,Ny,Nz);
  230. double *Ec = new double [Nx*Ny*Nz];
  231. InitEc(density,rs,Ec,Nx,Ny,Nz);
  232. cout<<"Correlation energy: "<<Integrate3(Ec,Nx,Ny,Nz)<<" Hartree"<<endl;
  233. double *left_Ec_m = new double [Nx*Ny*Nz];
  234. InitLeftEc_m(density,rs,left_Ec_m,Nx,Ny,Nz);
  235. double *right_Ec_m = new double [Nx*Ny*Nz];
  236. InitRightEc_m(density,rs,right_Ec_m,Nx,Ny,Nz);
  237. double *EXCwrtM = new double [Nx*Ny*Nz];
  238. FormIntM(homo2,lumo2,density,Ex,left_Ex_m,right_Ex_m,
  239.                              Ec,left_Ec_m,right_Ec_m,EXCwrtM,Nx,Ny,Nz);
  240. cout<<"Integral of EXC wrt M:"
  241.     <<Integrate3(EXCwrtM,Nx,Ny,Nz)<<" Hartree"<<endl;
  242. delete [] left_Ex_m;
  243. delete [] right_Ex_m;
  244. delete [] left_Ec_m;
  245. delete [] right_Ec_m;
  246. delete [] EXCwrtM;
  247. double *left_Ex_r = new double [Nx*Ny*Nz];
  248. InitEx(left_dens,rsleft,left_Ex_r,Nx,Ny,Nz);
  249. double *right_Ex_r = new double [Nx*Ny*Nz];
  250. InitEx(right_dens,rsright,right_Ex_r,Nx,Ny,Nz);
  251. double *left_Ec_r = new double [Nx*Ny*Nz];
  252. InitEc(left_dens,rsleft,left_Ec_r,Nx,Ny,Nz);
  253. double *right_Ec_r = new double [Nx*Ny*Nz];
  254. InitEc(right_dens,rsright,right_Ec_r,Nx,Ny,Nz);
  255. double *EXCwrtR = new double [Nx*Ny*Nz];
  256. FormIntR(homo2,lumo2,delta_dens,Ex,left_Ex_r,right_Ex_r,
  257.                                 Ec,left_Ec_r,right_Ec_r,EXCwrtR, Nx,Ny,Nz);
  258. cout<<"Integral of EXC wrt R:"
  259.     <<Integrate3(EXCwrtR,Nx,Ny,Nz)<<" Hartree"<<endl;
  260. delete [] density;
  261. delete [] left_dens;
  262. delete [] right_dens;
  263. delete [] delta_dens;
  264. delete [] rs;
  265. delete [] rsleft;
  266. delete [] rsright;
  267. delete [] EXCwrtR;
  268. delete [] left_Ex_r;
  269. delete [] right_Ex_r;
  270. delete [] left_Ec_r;
  271. delete [] right_Ec_r;
  272. delete [] Ex;
  273. delete [] Ec;
  274. double *diff = new double [Nx*Ny*Nz];
  275. Substract(homo2,lumo2,diff,Nx,Ny,Nz);
  276. delete [] homo2;
  277. delete [] lumo2;
  278. double *vpot = new double [Nx*Ny*Nz];
  279. void ConstructV(double *charge, double *dest, int Nx, int Ny, int Nz);
  280. ConstructV(diff,vpot,Nx,Ny,Nz);
  281. double *integralc = new double [Nx*Ny*Nz];
  282. Mult(diff,vpot,integralc,Nx,Ny,Nz);
  283. delete [] diff;
  284. delete [] vpot;
  285. cout<<"Coulomb Integral="<<Integrate3(integralc,Nx,Ny,Nz)<<endl;
  286. delete [] integralc;
  287. return 0;
  288. }/* END OF main() */
  289. void SquareMO(double *source,double *dest, int Nx, int Ny, int Nz)
  290. {/* To find norm of MO, see below */
  291.     for(int i=0;i<Nx;i+=1){
  292.         for(int j=0;j<Ny;j+=1){
  293.             for(int k=0;k<Nz;k+=1){
  294.                 dest[i*Ny*Nz+j*Nz+k]=source[i*Ny*Nz+j*Nz+k]
  295.                                     *source[i*Ny*Nz+j*Nz+k];
  296.             }
  297.         }
  298.     }
  299. }
  300. void Substract(double *term1, double *term2, double *source,
  301. int Nx, int Ny, int Nz)
  302. {
  303.     for(int i=0;i<Nx;i++){
  304.         for(int j=0;j<Ny;j++){
  305.             for(int k=0;k<Nz;k++){
  306.                 source[i*Ny*Nz+j*Nz+k]=term1[i*Ny*Nz+j*Nz+k]
  307.                                       -term2[i*Ny*Nz+j*Nz+k];
  308.             }
  309.         }
  310.     }
  311. }
  312. void Mult(double *term1, double *term2, double *source,
  313. int Nx, int Ny, int Nz)
  314. {
  315.     for(int i=0;i<Nx;i++){
  316.         for(int j=0;j<Ny;j++){
  317.             for(int k=0;k<Nz;k++){
  318.                 source[i*Ny*Nz+j*Nz+k]=term1[i*Ny*Nz+j*Nz+k]
  319.                                       *term2[i*Ny*Nz+j*Nz+k];
  320.             }
  321.         }
  322.     }
  323. }
  324. void InitDeltaDens(double *source,double *dest, int Nx, int Ny, int Nz)
  325. {
  326.     for(int i=0;i<Nx;i+=1){
  327.         for(int j=0;j<Ny;j+=1){
  328.             for(int k=0;k<Nz;k+=1){
  329.                 dest[i*Ny*Nz+j*Nz+k]=source[i*Ny*Nz+j*Nz+k]/100.;
  330.             }
  331.         }
  332.     }
  333. }
  334. void InitLeftDens(double *source,double *dest, double *diff,
  335. int Nx, int Ny, int Nz)
  336. {
  337.     for(int i=0;i<Nx;i+=1){
  338.         for(int j=0;j<Ny;j+=1){
  339.             for(int k=0;k<Nz;k+=1){
  340.                 dest[i*Ny*Nz+j*Nz+k]=source[i*Ny*Nz+j*Nz+k]
  341.                                     -diff[i*Ny*Nz+j*Nz+k];
  342.             }
  343.         }
  344.     }
  345. }
  346. void InitRightDens(double *source,double *dest, double *diff,
  347. int Nx, int Ny, int Nz)
  348. {
  349.     for(int i=0;i<Nx;i+=1){
  350.         for(int j=0;j<Ny;j+=1){
  351.             for(int k=0;k<Nz;k+=1){
  352.                 dest[i*Ny*Nz+j*Nz+k]=source[i*Ny*Nz+j*Nz+k]
  353.                                     +diff[i*Ny*Nz+j*Nz+k];
  354.             }
  355.         }
  356.     }
  357. }
  358. void Renorm(double *func, int Nx, int Ny, int Nz, double norm=1)
  359. {/* Call for normalizing MO density to 1, must calculate norm before */
  360.     for(int i=0;i<Nx;i+=1){
  361.         for(int j=0;j<Ny;j+=1){
  362.             for(int k=0;k<Nz;k+=1){
  363.                 func[i*Ny*Nz+j*Nz+k]=func[i*Ny*Nz+j*Nz+k]/norm;
  364.             }  
  365.         }
  366.     }
  367. }
  368. void InitRS(double *source,double *dest, int Nx, int Ny, int Nz)
  369. {/* Constructs rs from density */
  370.     for(int i=0;i<Nx;i+=1){
  371.         for(int j=0;j<Ny;j+=1){
  372.             for(int k=0;k<Nz;k+=1){
  373.                 dest[i*Ny*Nz+j*Nz+k]=pow(0.75/PI/source[i*Ny*Nz+j*Nz+k],1./3.);
  374.             }
  375.         }
  376.     }
  377. }
  378. void InitEx(double *rho, double *rs, double *dest, int Nx, int Ny, int Nz)
  379. {/* Exchange energy on grid */
  380.     for(int i=0;i<Nx;i+=1){
  381.         for(int j=0;j<Ny;j+=1){
  382.             for(int k=0;k<Nz;k+=1){
  383.                 dest[i*Ny*Nz+j*Nz+k]=AEX*rho[i*Ny*Nz+j*Nz+k]/rs[i*Ny*Nz+j*Nz+k];
  384.             }
  385.         }
  386.     }
  387. }
  388. void InitLeftEx_m(double *rho, double *rs, double *dest, int Nx, int Ny, int Nz)
  389. {/* Left part of Ex deriv wrt m */
  390.     for(int i=0;i<Nx;i+=1){
  391.         for(int j=0;j<Ny;j+=1){
  392.             for(int k=0;k<Nz;k+=1){
  393.                 dest[i*Ny*Nz+j*Nz+k]=AEX*0.5*(pow(1-h,4/3.)+pow(1+h,4/3.))
  394.                                     *rho[i*Ny*Nz+j*Nz+k]/rs[i*Ny*Nz+j*Nz+k];
  395.             }  
  396.         }
  397.     }
  398. }
  399. void InitRightEx_m(double *rho, double *rs, double *dest,
  400.                    int Nx, int Ny, int Nz)
  401. {/* Right part of Ex deriv wrt m */
  402.     for(int i=0;i<Nx;i+=1){
  403.         for(int j=0;j<Ny;j+=1){
  404.             for(int k=0;k<Nz;k+=1){
  405.                 dest[i*Ny*Nz+j*Nz+k]=AEX*0.5*(pow(1+h,4/3.)+pow(1-h,4/3.))
  406.                                     *rho[i*Ny*Nz+j*Nz+k]/rs[i*Ny*Nz+j*Nz+k];
  407.             }
  408.         }
  409.     }
  410. }
  411. void InitEc(double *rho, double *rs, double *dest, int Nx, int Ny, int Nz)
  412. {/* Correlation energy on grid */
  413.     for(int i=0;i<Nx;i+=1){
  414.         for(int j=0;j<Ny;j+=1){
  415.             for(int k=0;k<Nz;k+=1){
  416.                 dest[i*Ny*Nz+j*Nz+k]=(
  417.                 -2.*A_0*(1.+alpha1_0*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_0*
  418.                 (beta1_0*sqrt(rs[i*Ny*Nz+j*Nz+k])
  419.                 +beta2_0*rs[i*Ny*Nz+j*Nz+k]
  420.                 +beta3_0*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  421.                 +beta4_0*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  422.                 )))+
  423.                 -2.*A_1*(1.+alpha1_1*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_1*
  424.                 (beta1_1*sqrt(rs[i*Ny*Nz+j*Nz+k])
  425.                 +beta2_1*rs[i*Ny*Nz+j*Nz+k]
  426.                 +beta3_1*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  427.                 +beta4_1*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  428.                 )))*0+
  429.                 -2.*A_a*(1.+alpha1_a*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_a*
  430.                 (beta1_a*sqrt(rs[i*Ny*Nz+j*Nz+k])
  431.                 +beta2_a*rs[i*Ny*Nz+j*Nz+k]
  432.                 +beta3_a*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  433.                 +beta4_a*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  434.                 )))*0
  435.                 )*rho[i*Ny*Nz+j*Nz+k];
  436.             }
  437.         }
  438.     }
  439. }
  440. void InitLeftEc_m(double *rho, double *rs, double *dest,
  441.                   int Nx, int Ny, int Nz)
  442. {/* Left part of deriv wrt m */
  443.     for(int i=0;i<Nx;i+=1){
  444.         for(int j=0;j<Ny;j+=1){
  445.             for(int k=0;k<Nz;k+=1){
  446.                 dest[i*Ny*Nz+j*Nz+k]=(
  447.                 -2.*A_0*(1.+alpha1_0*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_0*
  448.                 (beta1_0*sqrt(rs[i*Ny*Nz+j*Nz+k])
  449.                 +beta2_0*rs[i*Ny*Nz+j*Nz+k]
  450.                 +beta3_0*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  451.                 +beta4_0*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  452.                 )))*(1-f_lz*pow(0-h,4))+
  453.                 -2.*A_1*(1.+alpha1_1*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_1*
  454.                 (beta1_1*sqrt(rs[i*Ny*Nz+j*Nz+k])
  455.                 +beta2_1*rs[i*Ny*Nz+j*Nz+k]
  456.                 +beta3_1*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  457.                 +beta4_1*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  458.                 )))*f_lz*pow(0-h,4)+
  459.                 -2.*A_a*(1.+alpha1_a*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_a*
  460.                 (beta1_a*sqrt(rs[i*Ny*Nz+j*Nz+k])
  461.                 +beta2_a*rs[i*Ny*Nz+j*Nz+k]
  462.                 +beta3_a*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  463.                 +beta4_a*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  464.                 )))*f_lz*(1-pow(0-h,4))/f2z)
  465.                 *rho[i*Ny*Nz+j*Nz+k];
  466.             }
  467.         }
  468.     }
  469. }
  470. void InitRightEc_m(double *rho, double *rs, double *dest,
  471.                    int Nx, int Ny, int Nz)
  472. {/* Right part of deriv wrt m */
  473.     for(int i=0;i<Nx;i+=1){
  474.         for(int j=0;j<Ny;j+=1){
  475.             for(int k=0;k<Nz;k+=1){
  476.                 dest[i*Ny*Nz+j*Nz+k]=(
  477.                 -2.*A_0*(1.+alpha1_0*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_0*
  478.                 (beta1_0*sqrt(rs[i*Ny*Nz+j*Nz+k])
  479.                 +beta2_0*rs[i*Ny*Nz+j*Nz+k]
  480.                 +beta3_0*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  481.                 +beta4_0*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  482.                 )))*(1-f_rz*pow(0+h,4))+
  483.                 -2.*A_1*(1.+alpha1_1*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_1*
  484.                 (beta1_1*sqrt(rs[i*Ny*Nz+j*Nz+k])
  485.                 +beta2_1*rs[i*Ny*Nz+j*Nz+k]
  486.                 +beta3_1*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  487.                 +beta4_1*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  488.                 )))*f_rz*pow(0+h,4)+
  489.                 -2.*A_a*(1.+alpha1_a*rs[i*Ny*Nz+j*Nz+k])*log(1.+1./(2*A_a*
  490.                 (beta1_a*sqrt(rs[i*Ny*Nz+j*Nz+k])
  491.                 +beta2_a*rs[i*Ny*Nz+j*Nz+k]
  492.                 +beta3_a*pow(rs[i*Ny*Nz+j*Nz+k],1.5)
  493.                 +beta4_a*pow(rs[i*Ny*Nz+j*Nz+k],p+1.)
  494.                 )))*f_rz*(1-pow(0+h,4))/f2z)
  495.                 *rho[i*Ny*Nz+j*Nz+k];
  496.             }
  497.         }
  498.     }
  499. }
  500. void FormIntM(double *rhoHOMO, double *rhoLUMO, double *rho,
  501.               double *Ex, double *Exleft, double *Exright,
  502.               double *Ec, double *Ecleft, double *Ecright, double *dest,
  503.               int Nx, int Ny, int Nz)
  504.  
  505. {/* EXC wrt M on grid  */
  506.     for(int i=0;i<Nx;i+=1){
  507.         for(int j=0;j<Ny;j+=1){
  508.             for(int k=0;k<Nz;k+=1){
  509.                 dest[i*Ny*Nz+j*Nz+k]=
  510.                 (rhoHOMO[i*Ny*Nz+j*Nz+k]+rhoLUMO[i*Ny*Nz+j*Nz+k])
  511.                 *(rhoHOMO[i*Ny*Nz+j*Nz+k]+rhoLUMO[i*Ny*Nz+j*Nz+k])
  512.                 *(Exleft[i*Ny*Nz+j*Nz+k]+Ecleft[i*Ny*Nz+j*Nz+k]
  513.                 -2*Ex[i*Ny*Nz+j*Nz+k]-2*Ec[i*Ny*Nz+j*Nz+k]
  514.                 +Exright[i*Ny*Nz+j*Nz+k]+Ecright[i*Ny*Nz+j*Nz+k])
  515.                 /rho[i*Ny*Nz+j*Nz+k]/rho[i*Ny*Nz+j*Nz+k]/h/h;
  516.             }
  517.         }
  518.     }          
  519. }
  520.  
  521. void FormIntR(double *rhoHOMO, double *rhoLUMO, double *drho,
  522.               double *Ex, double *Exleft, double *Exright,
  523.               double *Ec, double *Ecleft, double *Ecright,double *dest,
  524.               int Nx, int Ny, int Nz)
  525. {/* EXC wrt rho on grid */
  526.     for(int i=0;i<Nx;i+=1){
  527.         for(int j=0;j<Ny;j+=1){
  528.             for(int k=0;k<Nz;k+=1){
  529.                 dest[i*Ny*Nz+j*Nz+k]=
  530.                 (rhoHOMO[i*Ny*Nz+j*Nz+k]-rhoLUMO[i*Ny*Nz+j*Nz+k])
  531.                 *(rhoHOMO[i*Ny*Nz+j*Nz+k]-rhoLUMO[i*Ny*Nz+j*Nz+k])
  532.                 *(Exleft[i*Ny*Nz+j*Nz+k]+Ecleft[i*Ny*Nz+j*Nz+k]
  533.                 -2*Ex[i*Ny*Nz+j*Nz+k]-2*Ec[i*Ny*Nz+j*Nz+k]
  534.                 +Exright[i*Ny*Nz+j*Nz+k]+Ecright[i*Ny*Nz+j*Nz+k])
  535.                 /drho[i*Ny*Nz+j*Nz+k]/drho[i*Ny*Nz+j*Nz+k];
  536.             }
  537.         }
  538.     }
  539. }
  540. void WriteGCube(string filename, double *data,int Nx, int Ny, int Nz)
  541. {
  542.     ofstream outfile;
  543.     outfile.open(filename.c_str());
  544.     if (!outfile.is_open())
  545.     {
  546.         cout<<"Cannot open file "<<filename<<"."<<endl;
  547.         return;
  548.     }
  549.     outfile.precision(6);
  550.     outfile<<"User defined data"<<endl;
  551.     outfile<<"Comment line"<<endl;
  552.     outfile<<Nat<<" "<<xmin<<" "<<ymin<<" "<<zmin<<endl;
  553.     outfile<<Nx<<" "<<dx<<" 0.000000  0.000000"<<endl;
  554.     outfile<<Ny<<" "<<"0.000000 "<<dy<<" 0.000000"<<endl;
  555.     outfile<<Nz<<" "<<"0.000000  0.000000 "<<dz<<endl;
  556.     for(int i=0;i<Nat;i++)
  557.         outfile<<"  1   1.000000   0.000000   0.000000   0.000000"<<endl;
  558.     for (int i=0;i<Nx;i++){
  559.         for (int j=0;j<Ny;j++){
  560.             for (int k=0;k<Nz;k++){
  561.                 outfile<<"  "<<data[i*Ny*Nz+j*Nz+k]<<"  ";
  562.             if (k % 6 == 5)
  563.                 outfile<<endl;
  564.             }
  565.         outfile<<endl;
  566.         }
  567.     }
  568.     outfile.close();
  569. }
  570. void ConstructV(double *charge, double *dest, int Nx, int Ny, int Nz)
  571. {
  572.     double *temp = new double [Nx*Ny*Nz];
  573.     double fac=0.346873950;
  574.     for(int i=0;i<Nx;i++){
  575.         cout<<"Iteration "<<i<<"/"<<Nx-1<<"... "<<endl;
  576.         for(int j=0;j<Ny;j++){
  577.             for(int k=0;k<Nz;k++){
  578. /* Calculate V only on grid points which are used by Integrate3 */
  579.                 if ( ((i%2 == 1) && (j%2 == 0) && (k%2 == 0)) ||
  580.                      ((i%2 == 0) && (j%2 == 1) && (k%2 == 0)) ||
  581.                      ((i%2 == 0) && (j%2 == 0) && (k%2 == 1)) )
  582.                     {
  583.                         for(int ii=0;ii<Nx;ii+=1){
  584.                             for(int jj=0;jj<Ny;jj+=1){
  585.                                 for(int kk=0;kk<Nz;kk+=1){
  586.                                     if (i!=ii || j!=jj || k!=kk)
  587.                                         temp[ii*Ny*Nz+jj*Nz+kk]=
  588.                                         charge[ii*Ny*Nz+jj*Nz+kk]
  589.                                         /sqrt((i-ii)*(i-ii)
  590.                                              +(j-jj)*(j-jj)
  591.                                              +(k-kk)*(k-kk))/dx;
  592.                                     else
  593.                                         temp[ii*Ny*Nz+jj*Nz+kk]=
  594.                                         charge[ii*Ny*Nz+jj*Nz+kk]/dx/fac;
  595.                                 }
  596.                             }
  597.                         }
  598.             dest[i*Ny*Nz+j*Nz+k]=Integrate3(temp,Nx,Ny,Nz);}
  599.             else dest[i*Ny*Nz+j*Nz+k]=0;  
  600.             }
  601.         }
  602.     }
  603.     delete [] temp;
  604. }
  605. double Integrate3(double *func, int Nx, int Ny, int Nz)
  606. {/* Integrates on grid */
  607.     double result=0;
  608.     for(int i=0;i<Nx-1;i+=2){
  609.         for(int j=0;j<Ny-1;j+=2){
  610.             for(int k=0;k<Nz-1;k+=2){
  611.                 result+=func[(i+1)*Ny*Nz+j*Nz+k]/3.
  612.                        +func[i*Ny*Nz+(j+1)*Nz+k]/3.
  613.                        +func[i*Ny*Nz+j*Nz+(k+1)]/3.;
  614.             }
  615.         }
  616.     }
  617.     return result*8*dx*dy*dz;
  618. }
  619.  
  620. double integrate6(double *func, int Nx, int Ny, int Nz)
  621. {/* Not used right now, more expensive */
  622.     double result=0;
  623.     for(int i=1;i<Nx-1;i+=2){
  624.         for(int j=1;j<Ny-1;j+=2){
  625.             for(int k=1;k<Nz-1;k+=2){
  626.                 result+=func[(i+1)*Ny*Nz+j*Nz+k]/6.
  627.                       +func[(i-1)*Ny*Nz+j*Nz+k]/6.
  628.                       +func[i*Ny*Nz+(j+1)*Nz+k]/6.
  629.                       +func[i*Ny*Nz+(j-1)*Nz+k]/6.
  630.                       +func[i*Ny*Nz+j*Nz+(k+1)]/6.                        
  631.                       +func[i*Ny*Nz+j*Nz+(k-1)]/6.;
  632.             }
  633.         }
  634.     }
  635.     return result*8*dx*dy*dz;
  636. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!