Hubbard模型的计算程序

最新推荐文章于 2024-08-20 20:41:22 发布

diradical

最新推荐文章于 2024-08-20 20:41:22 发布

阅读量2.3k

点赞数

CC 4.0 BY-SA版权

分类专栏：理论化学计算文章标签： Hubbard 计算化学 C++

本文链接：https://blog.youkuaiyun.com/luozhen07/article/details/46000137

这篇博客介绍了如何利用Intel Math Kernel Library (MKL) 中的LAPACKE_dsyevd函数来求解Hubbard模型的本征能量和本征函数。编译过程依赖于正确配置对intel MKL的引用。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

只求完工而不求代码效率，现在终于写完这个程序了。接下来就要做代码的优化和并行化了——事实上需要优化的内容简直太多。

最后求解本征能量与本征函数的时候使用了 intel Math Kernel Library 提供的函数，即 lapack 中的 LAPACKE_dsyevd 函数。编译的时候需要正确调用 intel MKL 才可以编译完成。

/* 
  luozhen, 2015-05-25 
  Prof. Cgliu's Group, itcc, NJU 
  intel MKL linked as: -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread -lm 
*/

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <string>
#include <cstring>
#include <time.h>
#include <cstdlib>
#include <vector>
#include <iomanip>
#include "mkl_lapacke.h"

using namespace std;

//debug 
//定义class的方式效率并不高 
class config_wave
{
public:
	int *alpha;
	int *beta;
};
//--end debug

//根据输入文件名确定输出文件名
string FileName(string name)
{
	int site = name.find_last_of('.');
	if (site == 0)
	{
		name = "output";
	}
	else if(site > 0)
	{
		name.erase(site, name.size() - site + 1);
	}
	return name;
}

//获取当前系统日期与时间 
string _local_time()
{
    char buf[64];
	time_t curtime = time(NULL);
    tm *ptm = localtime(&curtime);
    sprintf(buf, "%d/%02d/%02d %02d:%02d:%02d", ptm->tm_year+1900, ptm->tm_mon+1,
    	ptm->tm_mday, ptm->tm_hour, ptm->tm_min, ptm->tm_sec);
    string _curtime(buf);
    return _curtime;
}

//获取当前主机的名称
string HostName()
{
	string hostname;
	//对于 Linux 系统
	ifstream hostname_file("/etc/hostname", ifstream::in);
	hostname_file >> hostname;
	return hostname;
}

//计算组合数的函数，递归方法实现
//int类型时计算上限是C(32,16)
int combinatorial(int n, int k)
{
	if(n < 0 || k < 0)
	{
		cerr<<"Error when calculating the combinatorial number of ("
			<<n<<", "<<k<<")!"<<endl;
		exit(10);
	}
	if(k > n)
	{
		return 0;
	}
	else if(k > n/2)
	{
		k = n - k;
	}
	else if(k == 1)
	{
		return n;
	}
	else if(k == 0)
	{
		return 1;
	}
	else
	{
		return combinatorial(n-1, k-1) + combinatorial(n-1, k);
	}
}

//由索引值生成01序列
//因为更高效的bitset不能在运行时动态生成，所以只能退而求其次使用数组 
int* GetBits(int index, int atoms)
{
	int *bits;
	bits = new int [atoms];
	for(int i = 0; i != atoms; ++ i)
	{
		bits[i] = 0;
	}
	int n = atoms - 1;
	int k = atoms / 2 - 1;
	int flag = combinatorial(n, k);
	for(int i = 0; i != atoms; ++ i)
	{
		if(index > flag)
		{
			index -= flag;
			-- n;
		}
		else
		{
			bits[i] = 1;
			-- n;
			-- k;
		}
		if(k < 0)
		{
			break;
		}
		flag = combinatorial(n, k);
	}
	return bits;
}

//由组态生成索引，对单种自旋的电子有效 
int GetIndex(int *_bits, int atoms)
{
	vector<int> _site;
	_site.push_back(0);
	for(int ix = 0; ix !=