/*
DFT++ is a density functional package developed by the research group
of Professor Tomas Arias
Copyright 1996-2003 Sohrab Ismail-Beigi
This file is part of DFT++.
DFT++ is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
DFT++ is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with DFT++; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
Please see the file CREDITS for a list of authors.
For academic users, we request that publications using results obtained with
this software reference
"New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).
and, if using the wavelet basis, further reference
"Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).
and
"Robust ab initio calculation of condensed matter: transparent convergence through
semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).
For your convenience, preprints of the above articles may be obtained from
http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/
//
// Sohrab Ismail-Beigi August 2 1999
//
// This is the main DFT++ executable.
// It handles electronic minimizations, calculation of forces, and
// ionic dynamics. Ionic relaxations are not working yet.
// It has been cleaned up a bit and the structure should be clear.
// The basic idea is to do any MPI setup if needed, then
// to try to read the input file, and then to do any minimizations/force/
// dynamics as needed, and finally to dump the current state and exit.
//
/* $Id: dft-two.cpp,v 1.1.2.92 2003/05/29 18:54:16 ivan Exp $ */
#include "header.h"
#include "parallel.h"
#include <time.h>
// Default input filename
#define DEFAULT_INPUTFILE "dft.in"
// This report level reports all the basic messages we might
// be interested to the log file. More or less anal options
// can be found in header.h
#define REPORT_LEVEL DFT_BASIC_LOG
// Prints a banner
void print_banner(char **argv,char *input_filename);
// Parses the command line arguments to main(). See down below.
void parse_command_line(int argc, char **argv,
char *input_filename,
char *log_filename);
void optest(Everything &e);
void Etest(Everything &e);
// HACK ALERT!!!
// This variable stores the square magnitude of Ygrad globally, so that the
// wavelet invL routing can access it.
double global_abs2Yg;
/***************************************************************************
* *
* Main Program: *
* *
***************************************************************************/
int
main(int argc,char**argv)
{
// Do any system dependent setup if needed
System::GlobalInit(&argc,&argv);
#ifdef DFT_PROFILING
// Activate timers and counters
timerActivateAll();
counterActivateAll();
// Turn initialization timer on
timerOn(0); // Initialization timer on
#endif
// The variables needed to do the work. Contains inside everything
// we will need!
Everything everything;
/**********************************************************
** **
** Reading command line, setting up logfile. **
** **
**********************************************************/
// initialize input/output files to defaults.
char input_filename[DFT_FILENAME_LEN] = DEFAULT_INPUTFILE;
char log_filename[DFT_FILENAME_LEN] = "";
// Parse command line and get input and log filenames.
parse_command_line(argc,argv,input_filename,log_filename);
// Open the logfile and make the system logger point to it
Output logfile(REPORT_LEVEL,log_filename);
System::global_log = &logfile;
/**********************************************************
** Setup section: read and parser the input file, and **
** the preform the various setup routines to prepare **
** everything. **
**********************************************************/
// Try to open the input file
dft_text_FILE *inputfile = dft_text_fopen(input_filename,"r");
if (inputfile == NULL)
die("\n%s: can't read '%s'. Aborting.\n\n",argv[0], input_filename);
// Print a welcome banner with useful information
print_banner(argv,input_filename);
// expand all includes in the input file.
inputfile->expand_include();
// parse input file and do the setup
parse(inputfile,everything);
dft_text_fclose(inputfile);
// ipd: hack to check it works - move it in everything.setup
// and put the actual variable names in there
// System::setup_scratch(3, 50, sizeof(complex));
// System::setup_scratch(3, 25, sizeof(complex));
// System::setup_scratch(3, 30, sizeof(complex));
// initial setup
everything.setup();
#ifdef DFT_PROFILING
timerOff(0); // turn off initialization timer.
timerOn(1); // turn on total computation timer.
#endif
//optest(everything); exit(1);
// Calculate the local pseudopotential vector if not reading it from file
//
// NOTE: this routine is basis dependent, and so
// has to be supplied by the basis developers
//
// if doing band structure we supply the self consistent potential
// which would include Vloc from the ions.
// change this once starting from charge densities is included.
if(everything.cntrl.band_structure_flag != TRUE)
Vlocpot(everything);
// calculate electron independent energies just once
calc_core_ewald_pulay_energies(everything);
//TEST the symemtries
// symmetrize_n(everything.elecvars.Vlocpot, &everything.symm);
// everything.elecvars.Vlocpot.write("Vs");
// here we may try to
if(everything.cntrl.finite_diff_test_flag) {
finitedifftest(everything);
die("FINITE DIFF TEST DONE, EXITTING\n");
}
//calc_all_energies(everything);
minimize_elec(everything);
//now do the forces - no iondynamics loop yet
if ( everything.cntrl.calculate_forces_flag ) {
calc_UVCn_d(everything);
// Calculate and report the forces
calc_ionic_forces(everything);
everything.ioninfo.print_force(dft_global_log, everything);
}
// Setup the signal handlers: turned off for now...
// the signal setup is disable right now. only a skeleton is left.
// setup_signals(&elecinfo,&elecvars);
#ifdef DFT_PROFILING
timerOff(1); // turn off total computation timer.
#endif
// Write out final electronic state to disk
dft_log("\nDone! Dumping final variables:\n\n");
dft_log_flush();
dump_and_stamp_variables(everything);
#ifdef DFT_PROFILING
// Print a final timer and counter report (with some very large
// iteration number, which should be our last!
timer_counter_Report(9999);
#endif
// End any system dependent stuff
System::GlobalFinalize();
return 0;
}
//
// Print the a banner with time, date, and parallelization info (if any)
//
void
print_banner(char **argv,char *input_filename)
{
time_t timenow = time(0);
dft_log("\n******************************************************\n");
dft_log("Current date and time: %s\n",ctime(&timenow));
#ifdef DFT_MPI
dft_log("This is an MPI run: N_Procs = %d\n",System::Get_N_Procs());
dft_log("Setting NODE_ID dependent seed\n");
srand(System::Get_procID()*1234+5678);//include time in it later
#endif
#ifdef DFT_THREAD
dft_log("This is a thread run: N_threads = %d\n",System::Get_N_threads());
#endif
dft_log("%s: input file is '%s'\n",argv[0],input_filename);
}
//
// Print out help information (e.g. command line options)
//
void
help(char *name)
{
// print out some help information
dft_log("\nUsage: %s [options]\n",name);
dft_log("\n\tPerforms ab-initio density-functional calculations.\n\n");
dft_log("options:\n\n");
dft_log("\t-h\t\thelp (this output)\n");
dft_log("\t-i <filename>\tspecifies input file, default = 'dft.in'\n");
dft_log("\t-o <filename>\tspecifies output log file, default = stdout\n");
#ifdef DFT_THREAD
dft_log("\t-n <integer>\tspecifies the number of threads, default = 1\n");
#endif
dft_log("\t-t\t\tprints an example input file to stdout and exits.\n");
dft_log("\n");
}
//
// Parse command line. Command line options currently
// are -i, -o, -t, and -h.
//
void
parse_command_line(int argc, char **argv,
char *input_filename,
char *log_filename)
{
// Logs to stdout. We need this to report
// basic errors and problems on the command line.
// Make the system logging point to it for now.
Output helplog;
System::global_log = &helplog;
int k=0;
while(1)
{
// If we get to the end of the command line with no problems,
// we are done!
if (++k==argc)
return;
// produce a template input file.
else if (MATCH(argv[k],"-t"))
{
Everything tempeverything;
print_default_template(dft_global_log,tempeverything);
break;
}
// override default output file.
else if (MATCH(argv[k],"-o"))
{
if (++k<argc)
strcpy(log_filename,argv[k]);
else
{
dft_log("\nOption -o requires an argument.\n\n");
break;
}
}
// override default input file.
else if (MATCH(argv[k],"-i"))
{
if (++k<argc)
strcpy(input_filename,argv[k]);
else
{
dft_log("\nOption -i requires an argument.\n\n");
break;
}
}
#ifdef DFT_THREAD
// override the default number of threads.
else if (MATCH(argv[k],"-n"))
{
if (++k<argc)
System::N_threads=atoi(argv[k]);
else
{
dft_log("\nOption -n requires an integer argument.\n\n");
break;
}
}
#endif
// Help or known option
else
{
// If the option is known, complain
if (!MATCH(argv[k],"-h"))
dft_log("\nUnknown option: %s.\n",argv[k]);
break;
}
}
// If we got here, print help and die
help(argv[0]);
System::GlobalFinalize();
exit(0);
}
// ipd: what is the argument type? template this? noooooooooooooooooooo -csg
// typedef rulez!!!!! ivan is a silly boy
// ipd: yeah dude, it rulez, but it's k&r, wanna switch to fortran???
void
compare(GenericColumn &c1, GenericColumn &c2,
char name1[20]="c1", char name2[20]="c2")
{
dft_log("\nComparing: %s (%s , %s) ",
name1,c1.getrepresentation(),c1.getembedding());
dft_log("and %s (%s , %s)\n",
name2,c2.getrepresentation(),c2.getembedding());
double m1, m2, m;
m1=sqrt(abs(c1^c1));
m2=sqrt(abs(c2^c2));
{
GenericColumn d(c1);
d -= c2;
m=sqrt(abs(d^d));
}
dft_log("|%s|\t= %le\n", name1, m1);
dft_log("|%s|\t= %le\n", name2, m2);
dft_log("|%s - %s|\t= %le\n", name1, name2, m);
// this is nasty, maybe we only want to print it if m1 and m2 are not small?
// but its not important...
complex ipd = (c1^c2);
complex nipd = ipd/(1.e-100+m1*m2);
dft_log("|%s^%s|\t= %le + %le i\n[|%s^%s|]\t= %le + %le i\n\n",
name1, name2, ipd.x, ipd.y, name1, name2, nipd.x, nipd.y);
}
void
optest(Everything &e)
{
int i, j, q;
dft_log("accessing via ColumnBundle.cols[i].data\n");
for (q = 0; q < e.elecvars.nstates; q++)
for(i = 0; i < e.elecvars.states[q].Y.my_ncols; i++)
for(j = 0; j < e.elecvars.states[q].Y.col[i]->length; j++){
(void)e.elecvars.states[q].Y.col[i]->data.d[j].x;
(void)e.elecvars.states[q].Y.col[i]->data.d[j].y;
}
RealSpaceScalarFieldColumn &n = e.elecvars.n;
RealSpaceScalarFieldColumn &d = e.elecvars.d;
CoeffSpaceScalarFieldColumn m(d), mmm(d);
RealSpaceScalarFieldColumn mm(n);
real mag = 0.0;
/**************************************
* LJ1 = 0
**************************************/
dft_log("=============================\n");
dft_log("LJ1 = 0 Test\n");
dft_log("=============================\n");
n.setmodones();
dft_log("Calling J(n)\n"); dft_log_flush();
apply_J(n,m);
dft_log("Calling L(m)\n"); dft_log_flush();
apply_L(m,mmm);
compare(m,mmm, "n", "LJn");
/**************************************
* IJ = 1
**************************************/
dft_log("=============================\n");
dft_log("IJ = 1 Test\n");
dft_log("=============================\n");
n.randomize();
dft_log("Calling J(n)\n"); dft_log_flush();
apply_J(n,m);
dft_log("Calling I(m)\n"); dft_log_flush();
apply_I(m,mm);
compare(n, mm, "n", "IJn");
/**************************************
* JI = 1
**************************************/
dft_log("=============================\n");
dft_log("JI = 1 Test\n");
dft_log("=============================\n");
d.randomize();
dft_log("Calling I(d)\n"); dft_log_flush();
apply_I(d,mm);
dft_log("Calling J(mm)\n"); dft_log_flush();
apply_J(mm,mmm);
compare(d, mmm, "d", "JId");
/**************************************
* Idag = I^
**************************************/
dft_log("=============================\n");
dft_log("Idag = I^ Test\n");
dft_log("=============================\n");
d.randomize();
dft_log("d: %s %s\n", d.getrepresentation(), d.getembedding());
m.randomize();
dft_log("m: %s %s\n", m.getrepresentation(), m.getembedding());
dft_log("Calling I(m)\n"); dft_log_flush();
apply_I(m,m);// must be in real space
dft_log("Calling I(d)\n"); dft_log_flush();
apply_I(d,mm);
dft_log("Calling Idag(m)\n"); dft_log_flush();
apply_Idag(m,mmm);
{
complex diff=(m^mm)/(mmm^d)-1;
dft_log("Idag is dag of I to %le + %le i\n\n",diff.x, diff.y );
}
/**************************************
* Jdag = J^
**************************************/
dft_log("=============================\n");
dft_log("Jdag = J^ Test\n");
dft_log("=============================\n");
d.randomize();
dft_log("d: %s %s\n", d.getrepresentation(), d.getembedding());
m.randomize();
dft_log("m: %s %s\n", m.getrepresentation(), m.getembedding());
dft_log("Calling J(m)\n"); dft_log_flush();
apply_J(m,mmm);// must be in real space
dft_log("Calling Jdag(d)\n"); dft_log_flush();
apply_Jdag(d,mm);
{
complex diff=(m^mm)/(mmm^d)-1;
dft_log("Jdag is dag of J to %le + %le i\n\n",diff.x, diff.y );
}
/**************************************
* J1^OJ1 = Vol
**************************************/
dft_log("=============================\n");
dft_log("J1^OJ1 = Vol Test\n");
dft_log("=============================\n");
for(i = 0; i < n.length; i++){
n.data.d[i].x = 1;
n.data.d[i].y = 0;
}
dft_log("Calling J(n)\n"); dft_log_flush();
apply_J(n,m);
dft_log("Calling O(m)\n"); dft_log_flush();
apply_O(m,mmm);
compare(m,mmm, "J1", "OJ1");
/**************************************
* J1^OJ1 = Vol on ColumnBundles
**************************************/
dft_log("==================================\n");
dft_log("(J1)^OJ1 = Vol on ColumnBundle Test\n");
dft_log("==================================\n");
for(i = 0; i < n.length; i++){
n.data.d[i].x = 1;
n.data.d[i].y = 0;
}
mag = 0.0;
for(i = 0; i < n.length; i++)
mag += abs(n.data.d[i]);
dft_log("\n mag[ n ] = %lf\n",mag);
ColumnBundle &C = e.elecvars.states[0].C;
ColumnBundle &Y = e.elecvars.states[0].Y;
for(i = 0; i < C.my_ncols; i++){
dft_log("Calling J(n)\n");
apply_J(n, *(C.col[i]));
}
dft_log("Calling O(Y)\n"); dft_log_flush();
O(C,Y);
dft_log("Calling C^Y\n"); dft_log_flush();
Matrix COC = C^Y;
dft_log("C^OC : \n");
COC.print();
dft_log("\n");
/**************************************
* Orthonorm test (on ColumnBundles)
**************************************/
dft_log("==================================\n");
dft_log("Orthonorm test (on ColumnBundles)\n");
dft_log("==================================\n");
C.randomize();
COC=C^O(C);
dft_log("\nOverlaps before:\n");
COC.print();
COC.hermetian = TRUE;
Matrix Umh=Uminusonehalf(COC,e.elecvars.states[0].W,e.elecvars.states[0].mu);
Y=C*Umh;
dft_log("\nOverlaps after:\n");
(Y^O(Y)).print();
dft_log("\n");
/**************************************
* Test on B*M algebra: (B*B)**2=B*B^
**************************************/
dft_log("==================================\n");
dft_log("B*M algebra test (B*B^)**2=B*B^ [idenpotency]\n");
dft_log("==================================\n");
//get idenpotent projector w/o O:
COC = C^C;
COC.hermetian = TRUE;
Y = C*Uminusonehalf(COC,e.elecvars.states[0].W,e.elecvars.states[0].mu);
dft_log("Projector idenpotency check:\n");
(Y^Y).print();
int col_num_diff = 3; // whatever to make it different
ColumnBundle B(Y.col_length, Y.tot_ncols-col_num_diff, Y.col[0]->basis);
// just put junk there, we're testing B*M
B.randomize();
ColumnBundle B1(B);
B1 = Y*(Y^B);
// project again
B = Y*(Y^B1);
dft_log("\n(PB)^(PB):\n");
(B^B).print();
//now compare
B-=B1;
dft_log("\nB^(PP-P)^(PP-P)B:\n");
(B^B).printe();
dft_log("\n");
}
void
Etest(Everything &e)
{
// calc_UVCn(e.elecinfo, e.elecvars, e.lattice, e.symm);
// solve_poisson(e.elecvars.n,e.elecvars.d);
calc_UVCn_d(e);
//e.elecvars.n.read("n_new");
//e.elecvars.n_up.read("n_up_new");
//e.elecvars.n_dn.read("n_dn_new");
calc_all_energies(e);
//e.elecvars.n_up.write("n_up");
//e.elecvars.n_dn.write("n_dn");
e.energies.print(System::global_log);
}
syntax highlighted by Code2HTML, v. 0.9.1