在openmx/source里的最重要的文件openmx.c。
下面简要罗列源码的思路:
1.定义CPU计算核数、MD计算迭代次数、统计时间用的一些临时变量。接下来进行MPI初始化,为并行计算做准备,openmx3.8版本只有并行版本。接着统计时间准备,便于监视计算耗时。最后检查输入变量。
2.确定openmx并行化使用的线程数,使用MPI_Bcast函数将来源于序列号为Host_ID的进程的消息广播发送到组内所有进程。
MPI_Bcast(&po, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
MPI_BCAST(buffer,count,datatype,root,comm)
IN/OUT buffer 通信消息缓冲区的起始地址(可变)
IN count 通信消息缓冲区中的数据个数(整型)
IN datatype 通信消息缓冲区中的数据类型(句柄)
IN root 发送广播的根的序列号(整型)
IN comm 通信子(句柄)
int MPI_Bcast(void* buffer,int count,MPI_Datatype datatype,int root, MPI_Comm comm)
设置完MPI后,进入判断全局参数的过程比如:
if (strcmp(argv[1],"-show")==0){
Show_DFT_DATA(argv);
exit(1);
}
判断是否进行NEB计算
if (neb_check(argv)) neb(argc,argv);
分配计算时间,并显示问候信息。
CompTime = (double**)malloc(sizeof(double*)*numprocs);
for (i=0; i<numprocs; i++){
CompTime[i] = (double*)malloc(sizeof(double)*30);
for (j=0; j<30; j++) CompTime[i][j] = 0.0;
}
if (myid==Host_ID){
printf("\n*******************************************************\n");
printf("*******************************************************\n");
printf(" Welcome to OpenMX Ver. %s \n",Version_OpenMX);
printf(" Copyright (C), 2002-2014, T. Ozaki \n");
printf(" OpenMX comes with ABSOLUTELY NO WARRANTY. \n");
printf(" This is free software, and you are welcome to \n");
printf(" redistribute it under the constitution of the GNU-GPL.\n");
printf("*******************************************************\n");
printf("*******************************************************\n\n");
}
Init_List_YOUSO();
remake_headfile = 0;
ScaleSize = 1.2;
接着读入输入文件
init_alloc_first();
在函数init_alloc_first()中定义了一个数组alloc_first[],每一个位对应一些选项,例如
alloc_first[26] = 1;
/***********************************************
Index_Snd_Grid_B2C
Index_Rcv_Grid_B2C
跟着文件读入的内容赋值:
CompTime[myid][1] = readfile(argv);
初始化内存路径:
sprintf(fileMemory,"%s%s.memory%i",filepath,filename,myid);
PrintMemory(fileMemory,0,"init");
PrintMemory_Fix();
进行全局初始化:
init();
开始DFT计算,
检查-mltest2模式,检查-forcetest2模式,检查强制一致性,代码见附录。
原胞和内部坐标优化:
/*******************************************************
optimization of cell and internal coordinates
*******************************************************/
if (CellOpt_switch!=0) cellopt(argv, CompTime);
进入核心计算部分自洽场-密度泛函(SCF-DFT)计算。在分子动力学仿真MD和几何优化(Geometry Optimization)中,前一步产生的重启文件(restart files)在后一步中可以加速借助于外插机制的收敛过程。被使用的之前产生的分子动力学或者几何优化的数目可以用关键字scf.ExtCharg.History来控制。
/****************************************************
SCF-DFT calculations, MD and geometrical
optimization.
****************************************************/
MD_iter = 1;
Temp_MD_iter = 1;
do {
if (MD_switch==12)
CompTime[myid][2] += truncation(1,1); /* EvsLC */
else if (MD_cellopt_flag==1)
CompTime[myid][2] += truncation(1,1); /* cell optimization */
else
CompTime[myid][2] += truncation(MD_iter,1);
if (ML_flag==1 && myid==Host_ID) Get_VSZ(MD_iter);
if (Solver==4) {
TRAN_Calc_GridBound( mpi_comm_level1, atomnum, WhatSpecies, Spe_Atom_Cut1,
Ngrid1, Grid_Origin, Gxyz, tv, gtv, rgtv, Left_tv, Right_tv );
/* output: TRAN_region[], TRAN_grid_bound */
}
if (Solver!=4 || TRAN_SCF_skip==0){
CompTime[myid][3] += DFT(MD_iter,(MD_iter-1)%orbitalOpt_per_MDIter+1);
iterout(MD_iter+MD_Current_Iter,MD_TimeStep*(MD_iter+MD_Current_Iter-1),filepath,filename);
/* MD or geometry optimization */
if (ML_flag==0) CompTime[myid][4] += MD_pac(MD_iter,argv[1]);
}
MD_iter++;
Temp_MD_iter++;
} while(MD_Opt_OK==0 && (MD_iter+MD_Current_Iter)<=MD_IterNumber);
if ( TRAN_output_hks ) {
/* left is dummy */
TRAN_RestartFile(mpi_comm_level1, "write","left",filepath,TRAN_hksoutfilename);
}
上面的过程包含了常规的密度泛函计算,后面的代码包括一些后处理和非平衡格林函数的计算。
附录:
全代码openmx.c
/*****************************************************************************
Ver. 3.7 (15/May/2013)
OpenMX (Open source package for Material eXplorer) is a program package
for linear scaling density functional calculations of large-scale materials.
Almost time-consuming parts can be performed in O(N) operations where N is
the number of atoms (or basis orbitals). Thus, the program might be useful
for studies of large-scale materials.
The distribution of this program package follows the practice of
the GNU General Public Licence (GPL).
OpenMX is based on
* Local density and generalized gradient approximation (LDA, LSDA, GGA)
to the exchange-corellation term
* Norm-conserving pseudo potentials
* Variationally optimized pseudo atomic basis orbitals
* Solution of Poisson's equation using FFT
* Evaluation of two-center integrals using Fourier transformation
* Evaluation of three-center integrals using fixed real space grids
* Simple mixing, direct inversion in the interative subspace (DIIS),
and Guaranteed-reduction Pulay's methods for SCF calculations.
* Solution of the eigenvalue problem using O(N) methods
* ...
See also our website (http://www.openmx-square.org/)
for recent developments.
**************************************************************
Copyright
Taisuke Ozaki
Present (15/May/2013) official address
Japan Advanced Institute of Science and Technology (JAIST)
Asahidai 1-1, Nomi, Ishikawa 923-1292, Japan
e-mail: [email protected]
**************************************************************
*****************************************************************************/
/**********************************************************************
openmx.c:
openmx.c is the main routine of OpenMX.
Log of openmx.c:
5/Oct/2003 Released by T.Ozaki
***********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/* stat section */
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
/* end stat section */
#include "openmx_common.h"
#include "mpi.h"
#include <omp.h>
#include "tran_prototypes.h"
#include "tran_variables.h"
int main(int argc, char *argv[])
{
static int numprocs,myid;
static int MD_iter,i,j,po,ip;
static char fileMemory[YOUSO10];
double TStime,TEtime;
/* MPI initialize */
mpi_comm_level1 = MPI_COMM_WORLD;
MPI_COMM_WORLD1 = MPI_COMM_WORLD;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
NUMPROCS_MPI_COMM_WORLD = numprocs;
MYID_MPI_COMM_WORLD = myid;
Num_Procs = numprocs;
/* for measuring elapsed time */
dtime(&TStime);
/* check argv */
if (argc==1){
printf("\nCould not find an input file.\n\n");
MPI_Finalize();
exit(0);
}
/* initialize Runtest_flag */
Runtest_flag = 0;
/****************************************************
./openmx -nt #
specifies the number of threads in parallelization
by OpenMP
****************************************************/
openmp_threads_num = 1; /* default */
po = 0;
if (myid==Host_ID){
for (i=0; i<argc; i++){
if ( strcmp(argv[i],"-nt")==0 ){
po = 1;
ip = i;
}
}
}
MPI_Bcast(&po, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
MPI_Bcast(&ip, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
if ( (argc-1)<(ip+1) ){
if (myid==Host_ID){
printf("cannot find the number of threads\n");
}
MPI_Finalize();
exit(0);
}
if ( po==1 ){
openmp_threads_num = atoi(argv[ip+1]);
if (openmp_threads_num<=0){
if (myid==Host_ID){
printf("check the number of threads\n");
}
MPI_Finalize();
exit(0);
}
}
omp_set_num_threads(openmp_threads_num);
if (myid==Host_ID){
printf("\nThe number of threads in each node for OpenMP parallelization is %d.\n\n",openmp_threads_num);
}
/****************************************************
./openmx -show directory
showing PAOs and VPS used in input files stored in
"directory".
****************************************************/
if (strcmp(argv[1],"-show")==0){
Show_DFT_DATA(argv);
exit(1);
}
/****************************************************
./openmx -maketest
making of *.out files in order to check whether
OpenMX normally runs on many platforms or not.
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketest")==0){
Maketest("S",argc,argv);
exit(1);
}
/****************************************************
./openmx -runtest
check whether OpenMX normally runs on many platforms
or not by comparing the stored *.out and generated
*.out on your machine.
****************************************************/
if ( strcmp(argv[1],"-runtest")==0){
Runtest("S",argc,argv);
}
/****************************************************
./openmx -maketestL
making of *.out files in order to check whether
OpenMX normally runs for relatively large systems
on many platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestL")==0){
Maketest("L",argc,argv);
exit(1);
}
/****************************************************
./openmx -maketestL2
making of *.out files in order to check whether
OpenMX normally runs for large systems on many
platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestL2")==0){
Maketest("L2",argc,argv);
exit(1);
}
/****************************************************
./openmx -runtestL
check whether OpenMX normally runs for relatively
large systems on many platforms or not by comparing
the stored *.out and generated *.out on your machine.
****************************************************/
if (strcmp(argv[1],"-runtestL")==0){
Runtest("L",argc,argv);
}
/****************************************************
./openmx -runtestL2
check whether OpenMX normally runs for large systems
on many platforms or not by comparing the stored *.out
and generated *.out on your machine.
****************************************************/
if (strcmp(argv[1],"-runtestL2")==0){
Runtest("L2",argc,argv);
}
/*******************************************************
check memory leak by monitoring actual used memory size
*******************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-mltest")==0){
Memory_Leak_test(argc,argv);
exit(1);
}
/****************************************************
./openmx -maketestG
making of *.out files in order to check whether
OpenMX normally runs for geometry optimization
on many platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestG")==0){
Maketest("G",argc,argv);
exit(1);
}
/****************************************************
./openmx -maketestC
making of *.out files in order to check whether
OpenMX normally runs for geometry optimization
on many platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestC")==0){
Maketest("C",argc,argv);
exit(1);
}
/****************************************************
./openmx -runtestG
check whether OpenMX normally runs for geometry
optimization on many platforms or not by comparing
the stored *.out and generated *.out on your machine.
****************************************************/
if (strcmp(argv[1],"-runtestG")==0){
Runtest("G",argc,argv);
}
/****************************************************
./openmx -runtestC
check whether OpenMX normally runs for simultaneous
optimization for cell and geometry on many platforms
or not by comparing the stored *.out and generated
*.out on your machine.
****************************************************/
if (strcmp(argv[1],"-runtestC")==0){
Runtest("C",argc,argv);
}
/****************************************************
./openmx -maketestWF
making of *.out files in order to check whether
OpenMX normally runs for generation of Wannier
functions on many platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestWF")==0){
Maketest("WF",argc,argv);
exit(1);
}
/****************************************************
./openmx -runtestWF
check whether OpenMX normally runs for generating
Wannier functions on many platforms or not by comparing
the stored *.out and generated *.out on your machine.
****************************************************/
if (strcmp(argv[1],"-runtestWF")==0){
Runtest("WF",argc,argv);
}
/****************************************************
./openmx -maketestNEGF
making of *.out files in order to check whether
OpenMX normally runs for NEGF calculations
on many platforms or not
****************************************************/
if ( (argc==2 || argc==3) && strcmp(argv[1],"-maketestNEGF")==0){
Maketest("NEGF",argc,argv);
exit(1);
}
/****************************************************
./openmx -runtestNEGF
check whether OpenMX normally runs for NEGF calculations
on many platforms or not
****************************************************/
if (strcmp(argv[1],"-runtestNEGF")==0){
Runtest("NEGF",argc,argv);
MPI_Finalize();
exit(1);
}
/*******************************************************
check consistency between analytic and numerical forces
*******************************************************/
if ( (argc==3 || argc==4) && strcmp(argv[1],"-forcetest")==0){
if (strcmp(argv[2],"0")==0) force_flag = 0;
else if (strcmp(argv[2],"1")==0) force_flag = 1;
else if (strcmp(argv[2],"2")==0) force_flag = 2;
else if (strcmp(argv[2],"3")==0) force_flag = 3;
else if (strcmp(argv[2],"4")==0) force_flag = 4;
else if (strcmp(argv[2],"5")==0) force_flag = 5;
else if (strcmp(argv[2],"6")==0) force_flag = 6;
else if (strcmp(argv[2],"7")==0) force_flag = 7;
else if (strcmp(argv[2],"8")==0) force_flag = 8;
else {
printf("unsupported flag for -forcetest\n");
exit(1);
}
Force_test(argc,argv);
exit(1);
}
/*******************************************************
check consistency between analytic and numerical stress
*******************************************************/
if ( (argc==3 || argc==4) && strcmp(argv[1],"-stresstest")==0){
if (strcmp(argv[2],"0")==0) stress_flag = 0;
else if (strcmp(argv[2],"1")==0) stress_flag = 1;
else if (strcmp(argv[2],"2")==0) stress_flag = 2;
else if (strcmp(argv[2],"3")==0) stress_flag = 3;
else if (strcmp(argv[2],"4")==0) stress_flag = 4;
else if (strcmp(argv[2],"5")==0) stress_flag = 5;
else if (strcmp(argv[2],"6")==0) stress_flag = 6;
else if (strcmp(argv[2],"7")==0) stress_flag = 7;
else if (strcmp(argv[2],"8")==0) stress_flag = 8;
else {
printf("unsupported flag for -stresstest\n");
exit(1);
}
Stress_test(argc,argv);
MPI_Finalize();
exit(1);
}
/*******************************************************
check the NEB calculation or not, and if yes, go to
the NEB calculation.
*******************************************************/
if (neb_check(argv)) neb(argc,argv);
/*******************************************************
allocation of CompTime and show the greeting message
*******************************************************/
CompTime = (double**)malloc(sizeof(double*)*numprocs);
for (i=0; i<numprocs; i++){
CompTime[i] = (double*)malloc(sizeof(double)*30);
for (j=0; j<30; j++) CompTime[i][j] = 0.0;
}
if (myid==Host_ID){
printf("\n*******************************************************\n");
printf("*******************************************************\n");
printf(" Welcome to OpenMX Ver. %s \n",Version_OpenMX);
printf(" Copyright (C), 2002-2014, T. Ozaki \n");
printf(" OpenMX comes with ABSOLUTELY NO WARRANTY. \n");
printf(" This is free software, and you are welcome to \n");
printf(" redistribute it under the constitution of the GNU-GPL.\n");
printf("*******************************************************\n");
printf("*******************************************************\n\n");
}
Init_List_YOUSO();
remake_headfile = 0;
ScaleSize = 1.2;
/****************************************************
Read the input file
****************************************************/
init_alloc_first();
CompTime[myid][1] = readfile(argv);
MPI_Barrier(MPI_COMM_WORLD1);
/* initialize PrintMemory routine */
sprintf(fileMemory,"%s%s.memory%i",filepath,filename,myid);
PrintMemory(fileMemory,0,"init");
PrintMemory_Fix();
/* initialize */
init();
/* for DFTD-vdW by okuno */
/* for version_dftD by Ellner*/
if(dftD_switch==1 && version_dftD==2) DFTDvdW_init();
if(dftD_switch==1 && version_dftD==3) DFTD3vdW_init();
/* check "-mltest2" mode */
po = 0;
if (myid==Host_ID){
for (i=0; i<argc; i++){
if ( strcmp(argv[i],"-mltest2")==0 ){
po = 1;
ip = i;
}
}
}
MPI_Bcast(&po, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
MPI_Bcast(&ip, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
if ( po==1 ) ML_flag = 1;
else ML_flag = 0;
/* check "-forcetest2" mode */
po = 0;
if (myid==Host_ID){
for (i=0; i<argc; i++){
if ( strcmp(argv[i],"-forcetest2")==0 ){
po = 1;
ip = i;
}
}
}
MPI_Bcast(&po, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
MPI_Bcast(&ip, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
if ( po==1 ){
force_flag = atoi(argv[ip+1]);
ForceConsistency_flag = 1;
}
/* check force consistency
the number of processes
should be less than 2.
*/
if (ForceConsistency_flag==1){
Check_Force(argv);
CompTime[myid][20] = OutData(argv[1]);
Merge_LogFile(argv[1]);
Free_Arrays(0);
MPI_Finalize();
exit(0);
return 1;
}
/* check "-stresstest2" mode */
po = 0;
if (myid==Host_ID){
for (i=0; i<argc; i++){
if ( strcmp(argv[i],"-stresstest2")==0 ){
po = 1;
ip = i;
}
}
}
MPI_Bcast(&po, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
MPI_Bcast(&ip, 1, MPI_INT, Host_ID, MPI_COMM_WORLD1);
if ( po==1 ){
stress_flag = atoi(argv[ip+1]);
StressConsistency_flag = 1;
}
/* check stress consistency
the number of processes
should be less than 2.
*/
if (StressConsistency_flag==1){
Check_Stress(argv);
CompTime[myid][20] = OutData(argv[1]);
Merge_LogFile(argv[1]);
Free_Arrays(0);
MPI_Finalize();
exit(0);
return 1;
}
/*******************************************************
optimization of cell and internal coordinates
*******************************************************/
if (CellOpt_switch!=0) cellopt(argv, CompTime);
/****************************************************
SCF-DFT calculations, MD and geometrical
optimization.
****************************************************/
MD_iter = 1;
Temp_MD_iter = 1;
do {
if (MD_switch==12)
CompTime[myid][2] += truncation(1,1); /* EvsLC */
else if (MD_cellopt_flag==1)
CompTime[myid][2] += truncation(1,1); /* cell optimization */
else
CompTime[myid][2] += truncation(MD_iter,1);
if (ML_flag==1 && myid==Host_ID) Get_VSZ(MD_iter);
if (Solver==4) {
TRAN_Calc_GridBound( mpi_comm_level1, atomnum, WhatSpecies, Spe_Atom_Cut1,
Ngrid1, Grid_Origin, Gxyz, tv, gtv, rgtv, Left_tv, Right_tv );
/* output: TRAN_region[], TRAN_grid_bound */
}
if (Solver!=4 || TRAN_SCF_skip==0){
CompTime[myid][3] += DFT(MD_iter,(MD_iter-1)%orbitalOpt_per_MDIter+1);
iterout(MD_iter+MD_Current_Iter,MD_TimeStep*(MD_iter+MD_Current_Iter-1),filepath,filename);
/* MD or geometry optimization */
if (ML_flag==0) CompTime[myid][4] += MD_pac(MD_iter,argv[1]);
}
MD_iter++;
Temp_MD_iter++;
} while(MD_Opt_OK==0 && (MD_iter+MD_Current_Iter)<=MD_IterNumber);
if ( TRAN_output_hks ) {
/* left is dummy */
TRAN_RestartFile(mpi_comm_level1, "write","left",filepath,TRAN_hksoutfilename);
}
/****************************************************
calculate Voronoi charge
****************************************************/
if (Voronoi_Charge_flag==1) Voronoi_Charge();
/****************************************************
calculate Voronoi orbital magnetic moment
****************************************************/
if (Voronoi_OrbM_flag==1) Voronoi_Orbital_Moment();
/****************************************************
output analysis of decomposed energies
****************************************************/
if (Energy_Decomposition_flag==1) Output_Energy_Decomposition();
/****************************************************
making of a file *.frac for the fractional coordinates
****************************************************/
Make_FracCoord(argv[1]);
/****************************************************
generate Wannier functions added by Hongming Weng
****************************************************/
/* hmweng */
if(Wannier_Func_Calc){
if (myid==Host_ID) printf("Calling Generate_Wannier...\n");fflush(0);
Generate_Wannier(argv[1]);
}
/*********************************************************
Electronic transport calculations based on NEGF:
transmission, current, eigen channel analysis, and
real space analysis of current
*********************************************************/
if (Solver==4 && TRAN_analysis==1) {
/* if SCF is skipped, calculate values of basis functions on each grid */
if (1<=TRAN_SCF_skip) i = Set_Orbitals_Grid(0);
if (SpinP_switch==3) {
TRAN_Main_Analysis_NC( mpi_comm_level1, argc, argv, Matomnum, M2G,
GridN_Atom, GridListAtom, CellListAtom,
Orbs_Grid, TNumGrid );
}
else {
TRAN_Main_Analysis( mpi_comm_level1, argc, argv, Matomnum, M2G,
GridN_Atom, GridListAtom, CellListAtom,
Orbs_Grid, TNumGrid );
}
}
/****************************************************
Making of output files
****************************************************/
if (OutData_bin_flag)
CompTime[myid][20] = OutData_Binary(argv[1]);
else
CompTime[myid][20] = OutData(argv[1]);
/****************************************************
write connectivity, Hamiltonian, overlap, density
matrices, and etc. to a file, filename.scfout
****************************************************/
if (HS_fileout==1) SCF2File("write",argv[1]);
/* elapsed time */
dtime(&TEtime);
CompTime[myid][0] = TEtime - TStime;
Output_CompTime();
for (i=0; i<numprocs; i++){
free(CompTime[i]);
}
free(CompTime);
/* merge log files */
Merge_LogFile(argv[1]);
/* free arrays for NEGF */
if (Solver==4){
TRAN_Deallocate_Atoms();
TRAN_Deallocate_RestartFile("left");
TRAN_Deallocate_RestartFile("right");
}
/* free arrays */
Free_Arrays(0);
/* print memory */
PrintMemory("total",0,"sum");
MPI_Barrier(MPI_COMM_WORLD);
if (myid==Host_ID){
printf("\nThe calculation was normally finished.\n");
}
MPI_Finalize();
return 0;
}