共价模型单种自旋的电子组态的索引方法

       此方法来自武剑的论文，见 Jian Wu and Yuansheng Jiang, J. Comp. Chem. 2000, 21, 856.

       此处的main函数只是用于测试算法是否正确，此算法最后会作为另一个程序的子函数。

/*
  luozhen, 2015-05-20
*/

//Available only for Sz = 0

#include <iostream>
#include <vector>
#include <bitset>

using namespace std;

//define the total number of atoms
const int atoms = 6;

//combinatorial number calculation
//compiled to 64bit binary, the maxinum number under handle is C(32,16)
int _combinatorial(int n, int k)
{
	if(k > n)
		return 0;
	if(k > n/2)
		k = n - k;
	if(k == 1) {
		return n;
	} else if(k == 0) {
		return 1;
	} else {
		return _combinatorial(n-1, k-1) + _combinatorial(n-1, k);
	}
}

//Applied to Alpha electrons
int GetAlphaIndex(bitset<atoms> _wave)
{
	vector<int> alpha_site;
	alpha_site.push_back(0);
	for(int ix = 0; ix != atoms; ++ ix) {
		if(_wave[ix]) {
			alpha_site.push_back(ix+1);
		}
	}
	int index = 1;
	for(int i = 1; i != atoms/2+1; ++ i) {
		if(alpha_site[i] - alpha_site[i-1] > 1) {
			for(int j = 1; j != alpha_site[i] - alpha_site[i-1]; ++ j) {
				index += _combinatorial(atoms - alpha_site[i-1] -j, atoms/2 - i);
			}
		}
	}
	return index;
}

//to check whether the code is right
int main()
{
	bitset<atoms> iwave;
	iwave.reset();
	int site1, site2, site3;
	cout<<"Enter site1, site2 and site3 (1 <= site_i <= "<<atoms<<"):"<<endl;
	cin>>site1>>site2>>site3;
	iwave.set(site1-1);
	iwave.set(site2-1);
	iwave.set(site3-1);
	int temp = GetAlphaIndex(iwave);
	cout<<"config: "<<site1<<site2<<site3
		<<"\tbitset: "<<iwave<<"\tindex: "<<GetAlphaIndex(iwave)<<endl;
	return 0;
}