НЕКОТОРЫЕ МОДУЛИ, НАПИСАННЫЕ НА С++

Здесь очень кратко описаны две небольшие программы, каждая из которых решает какую-то одну проблему молекулярной спектроскопии.

1. Вывод аналитического выражения для суммы диагональных элементов матрицы эффективного вращательного гамильтониана

Одной из задач молекулярной спектроскопии является определение высокоточной структуры молекулы по ее колебательно-вращательному спектру, то есть в преобразовании информации о молекуле, заключенной в частотах линий поглощения, в набор длин связей и углов между ними и связывающих атомы сил. Это задача в настоящее время решена только для двухатомных и некотороых трехатомных молекул. Для молекул с большим числом атомов задача сильно усложняется, так как требует очень большого числа экспериментальных данных и больших объемов вычислений, длящихся месяцы, а иногда и годы. Задачу условно можно разбить на нескольких этапов. Один из них состоит в определении из измеренного спектра параметров эффективного вращательного гамильтониана (оператора энергии), с помошью которого строится теоретическая модель наблюдаемого спектра. Не вдаваясь в специфические подробности, понятные только для узкого круга специалистов, следует сказать, что параметры эффективного гамильтониана сильно коррелированы, так как гамильтониан определен с точностью до произвольного унитарного преобразования. Поэтому параметры определяются с очень низкой точностью. В то же время имеются инварианты гамильтониана, остающиеся неизменными при любых унитарных преобразованиях. Они являются функциями параметров, малокоррелированы, и могут быть определены из спектра с высокой точностью. Таким образом, необходимо знать функциональную связь между используемыми инвариантами и параметрами гамильтониана.

Для получения такой функциональной зависимости в аналитическом виде, то есть по-сути вывода формулы, была написана программа на C++, представленная ниже. Для пояснения действий, которые выполняет данная программа, прежде всего остановимся на самом эффективном гамильтониане (для простейших квазижестких молекул типа асимметричного волчка). С точностью до поправок второго порядка он имеет следующий вид:

Heff
(1)

Здесь Baa и tau_abgd параметры гамильтониана. Ja декартовые компоненты оператора углового момента.

Оператор Гамильтона может быть записан в дифференциальной форме, как в выражении (1), так и переведен в матричную форму. Инвариантом(ами) гамильтониана, о котором(ых) речь шла выше, является сумма диагональных элементов его матрицы. Обычно для записи матрицы гамильтониана в качестве базиса используются функции симметричного волчка psi_jkm. В этом базисе матрицы компонент оператора углового момента имеют диагональный вид по квантовому числу j. Ненулевые матричные элементы в символах бра кет выглядят следующим образом:

(2)
(3)

Видно, что j правой и левой функций в интегралах для ненулевых матричных элементов должно совпадать. Поэтому матрица гамльтониана приобретает блочно диагональный вид с блоками, характеризующимися конкретными значениями j. Отметим, что квантовое число j может принимать любые положительные целочисленные значения. Для каждого блока может быть вычислена своя сумма диагональных элементов, которая есть инвариантом.

Результатом работы программы по выводу аналитического выражения для суммы диагональных элементов любого блока гамильтониана (1) с матричными элементами оператора углового момента (2) и (3) является следующее. Сумма диагональных элементов диагональных блоков с определенным j равна

(4)

Вывод этого соотношения основывался на формуле произведения матриц. Так для квартичной части гамильтониана по операторам углового момента использовалась следующая формула, дающая сумму диагональных элементов отдельных членов:

,
(5)

где

,

В эту формулу подставлялись выражения (2) и (3) для отдельных матричных элементов операторов углового момента.

Программа выполнена как проект console application в VC++ 6. Она состоит из следующих исполняемых и заголовочных модулей:

  • types.h
  • functions.h
  • Jway.h
  • StdAfx.h
  • main.cpp
  • functions.cpp
  • Jway.cpp

types.h Определяет составные типы данных, используемые в данной задаче.

#if !defined(types_h)
#define types_h

#include <vector>
#include <string>

typedef enum {x,y,z} rlist;
typedef struct
{
	int dir[4];  // direction
	int km[4]; // shift in K (dk)
} path;

// Информация о всех диагональных полных путях для одного квартичного 
// члена JaJbJcJd.
// Каждый путь представляет собой произведение четырех элементарных
// путей.
class Pathes
{
public:
	int op_r[4];           // Тип множителей Ja.
	std::vector<path> onepath;  // Направление элементарных путей и K для 
						   // каждого элементарного пути.
};

//typedef struct
class j_k_powers
{
public:
	double val; // factor;
	int jp;  // power of j(j+1);
	int kp;  // power of k;
    bool operator==(const j_k_powers& ) const;
};
//} j_k_powers;

typedef  std::vector<j_k_powers> Term;
typedef  std::vector<Term> Tape;

typedef struct
{
	double v;  // multiplier.
	int pi;  // power of i (sqrt(-1)).
} valim;
typedef std::vector<valim> allvm;

// Class finres should collect final results;
class finres
{
public:
	Term d;
	std::string sq;
	int pi;
	int used;
};

#endif //types_h

functions.h Объявляет прототипы функций.

#include "types.h"
#include <fstream>

int binom (int , int );
int ntape(int,int);
void Fill_tapes(Tape&, Tape&);
void print_Tape(Tape ,char* );
void print_Term(Term ,char* );
Term Terms_multiplication(Term ,Term );
bool compj_k_powers (j_k_powers, j_k_powers);
void EqualAdds(Term& );
Term Terms_summation(Term ,Term );
char fs(int);

Jway.h Объявляет класс Jway. Этот класс содержит "направления" и количество матричных элементов для базовых операторов углового момента Jx, Jy, Jz.

// Jway.h: interface for the Jway class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_JWAY_H__6F14C1BE_B8F4_443F_939E_07E831AB85C4__INCLUDED_)
#define AFX_JWAY_H__6F14C1BE_B8F4_443F_939E_07E831AB85C4__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
#include <vector>
#include "types.h"

class Jway  
{
public:
	Jway();
	virtual ~Jway();
	void Init(rlist);
	int nways;               // Number of the elementary pathes (1 or 2).
	std::vector<int> way;      //They tell about dk for the quantum path.
};
#endif // !defined(AFX_JWAY_H__6F14C1BE_B8F4_443F_939E_07E831AB85C4__INCLUDED_)

main.cpp Модуль, реализующий основной алгоритм.

#pragma warning(disable: 4786)  // turn off the debugger warming C4786.
#include "stdafx.h"
#include "functions.h"
#include "Jway.h"
#include <fstream>
#include <iostream.h>
#include <iomanip.h>
#include <vector>
#include <string>

using namespace std;
ofstream fn;

// =========  main =========================================
int main(int argc, char* argv[])
{
	fn.open("res.txt",ios::out);
    
	Tape pjk_tape, pJz;
	
	int i1,i2,i3,i4;
	char sj[5]; sj[4]='\0';
	Jway J[3];
	path w;
	Pathes setpathes;
	Tape onepath(1);
	Term d(1);
	Tape allpathes;
	valim vi;
	allvm pvm;
	vector<finres> fr;
	fr.reserve(100);

	Fill_tapes(pjk_tape, pJz);

	print_Tape(pjk_tape, "Matrix elements of J+");
	print_Tape(pJz, "Matrix elements of Jz");

	J[x].Init(x);
	J[y].Init(y);
	J[z].Init(z);

// Determination of the terms having diagonal elements.
// Determination of the ways from left state to state on the right.

	fn << "Степень " << "Степень " << "Свободный" << endl;
	fn << " J**2 " << "    Jz   " << "  член" << endl;
	int a=0;
	int m=0;

// Начало перечисления всех операторов типа JaJbJcJd
	for(i1=x; i1 <= z; i1++)
		for(i2=x; i2 <= z; i2++)
			for(i3=x; i3 <= z; i3++)
				for(i4=x; i4 <= z; i4++)
				{
// Names of dekart components.
					sj[0]=fs(i1);
					sj[1]=fs(i2);
					sj[2]=fs(i3);
					sj[3]=fs(i4);
//
					setpathes.op_r[0]=i1;
					setpathes.op_r[1]=i2;
					setpathes.op_r[2]=i3;
					setpathes.op_r[3]=i4;
					int id1=0;
					int np=0;
					int y=0;
					setpathes.onepath.clear();
					allpathes.clear();
// Loop through all pathes for a given term JaJbJcJd.
// Данные о путях, соответствующих диагональным элементам, записываются в
// векторе setpathes (Направления элементарных путей и значение изменения
// квантового числа K по сравнению с исходным).
		for(int k1=0; k1<J[i1].nways; k1++)
		{
			id1=J[i1].way[k1];
			int ip1=ntape(id1,J[i1].way[k1]);
			int id2;
			for(int k2=0; k2<J[i2].nways; k2++)
			{
				id2=id1+J[i2].way[k2];
				int ip2=ntape(id2,J[i2].way[k2]);
				int id3;
				for(int k3=0; k3<J[i3].nways; k3++)
				{
					id3=id2+J[i3].way[k3];
					int ip3=ntape(id3,J[i3].way[k3]);
					int id4;
					for(int k4=0; k4<J[i4].nways; k4++)
					{
						id4=id3+J[i4].way[k4];
						int ip4=ntape(id4,J[i4].way[k4]);

// Total way must be of deltK=2=2

						if(id4==2)
						{
							m++;  np++;  y=1;
							w.dir[0]=J[i1].way[k1];
							w.dir[1]=J[i2].way[k2];
							w.dir[2]=J[i3].way[k3];
							w.dir[3]=J[i4].way[k4];
							w.km[0]=ip1;
							w.km[1]=ip2;
							w.km[2]=ip3;
							w.km[3]=ip4;
// В setpathes запоминается информация о всех deltK=2 путях для данного члена JaJbJcJd.
							setpathes.onepath.push_back(w);
						}
						a++;
					}
				}
			}
		}
		if(y==1)  // There is a deltK=2 matrix element.
		{
			pvm.clear();
// Loop through the all pathes found for a given term JaJbJcJd.
			for(int i=0; i<setpathes.onepath.size(); i++)
			{
// One path multipliers and imaginary;
// Multiplication of the coefficients of elementary pathes.
				double v=1;  // multiplier.
				int pi=0;  // power of i, i.e. sqrt(-1).
				for(k=0; k<4; k++)
				{
					switch(setpathes.op_r[k]) // What is it? x, y, z?
					{
					case(0):   // x
						v*=0.5;
						break;
					case(1):  // y
						pi++;
						v*=0.5;
						if(setpathes.onepath[i].dir[k]==-1)
							v=-v;
						break;
					default:
						break;
					}
				}
// One path imaginary;
				switch(pi)
				{
				case(0):
					break;
				case(1):     // i
					break;
				case(2):     // i^2
					pi=0;
					v=-v;
					break;
				case(3):     // i^3
					pi=1;
					v=-v;
					break;
				case(4):     // i^4
					pi=0;
					break;
				}
// At this point pi=0 - real value, pi=1 - imaginary value;
				vi.v=v;
				vi.pi=pi;
				pvm.push_back(vi);

// It should be two raised uncompensated elementary pathes.
// Other two ones contain or lowered and raised pathes or two z-pathes.
// Choosing the raised and horisontal elementary pathes for one way.
// Vector onepath will contain analitical expressions for 
// unsquared part of matrix elements.
				onepath.clear();
				for(k=0; k<4; k++)
				{
					switch(setpathes.onepath[i].dir[k])
					{
					case (0):		// Jz path
						d=pJz[setpathes.onepath[i].km[k]+2];
						onepath.push_back(d);
						break;
					case (-1):		// Jx, Jy - lowered path
						d=pjk_tape[setpathes.onepath[i].km[k]+1];
						onepath.push_back(d);
						break;
					}
				}
				print_Tape(onepath, "Raised pathes");
				d=onepath[0];
// multiplication of the elementary pathes.
// Here we derive an anlaitical expression for a single path.
				for(k=1; k<onepath.size(); k++)
					d=Terms_multiplication(d,onepath[k]);
				print_Term(d, "One path");
				EqualAdds(d);
				print_Term(d, "Dropped with 0 value");
				allpathes.push_back(d);
			}
			print_Tape(allpathes,"allpathes");
// Correction of the coefficients
			for(i=0; i<setpathes.onepath.size(); i++)
				for(int k=0; k<allpathes[i].size(); k++)
					allpathes[i][k].val*=pvm[i].v;
			print_Tape(allpathes,"allpathes corrected");
// Final summation
			d=allpathes[0];
			for(int k=1; k<allpathes.size(); k++)
				d=Terms_summation(d,allpathes[k]);
			print_Term(d,"Final matrix element");
			finres fa;
			fa.sq=sj;
			fa.pi=vi.pi;
			fa.d=d;
			fa.used=-1;
			fr.push_back(fa);
			fn << fa.sq << endl;
//			fn << fa.pi << endl;
			if(fa.pi==0)
				fn << fa.pi << " (Real value)" << endl;
			else
				fn << fa.pi << " (Imaginary value)" << endl;
			for(i=0; i<fa.d.size(); i++)
			{
				fn.width(4);
				fn << fa.d[i].jp;
				fn.width(4);
				fn << fa.d[i].kp;
				fn.width(10);
				fn << fa.d[i].val << endl;
			}
		}
	}

// Final summation
	int k=0;
	int n=0;
	while( k < fr.size())
	{
		if(fr[k].used == -1)
		{
			fr[k].used=n;
			d=fr[k].d;
			int k1=k;
			string sw=fr[k].sq;
// Looking for terms of the same type
			if(k1+1 == fr.size())break;
			for(int m=k1+1; m < fr.size(); m++)
			{
				if(fr[m].used != -1)continue;
				string st=fr[m].sq;
				if((sw[0] == st[1] && sw[1] == st[0] && sw[2] == st[2] && sw[3] == st[3]) ||
				   (sw[0] == st[1] && sw[1] == st[0] && sw[2] == st[3] && sw[3] == st[2]) ||
				   (sw[0] == st[0] && sw[1] == st[1] && sw[2] == st[3] && sw[3] == st[2]) ||
				   (sw[0] == st[2] && sw[1] == st[3] && sw[2] == st[0] && sw[3] == st[1]) ||
				   (sw[0] == st[2] && sw[1] == st[3] && sw[2] == st[1] && sw[3] == st[0]) ||
				   (sw[0] == st[3] && sw[1] == st[2] && sw[2] == st[0] && sw[3] == st[1]) ||
				   (sw[0] == st[3] && sw[1] == st[2] && sw[2] == st[1] && sw[3] == st[0]))
				{
					fr[m].used=n;
					d=Terms_summation(d,fr[m].d);
				}
			}
			print_Term(d,"");
			n++;
		}
		k++;
	}

    cout << a << endl;
    cout << m << endl;
	return 0;
}

functions.cpp Часть кода, вынесенного в функции.


#pragma warning(disable: 4786)  // turn off the debugger warming C4786.
#include "stdafx.h"
#include "functions.h"
#include <vector>
#include "types.h"
#include <algorithm>
#define wmn 1000

using namespace std;

extern ofstream ff;
// ntape определяет полосу, в которой находится элементарный путь.
// Если путь восходящий, то возвращается левое K матричного элемента.
// Если путь нисходящий, то возвращается правое K матричного элемента.
// Если путь горизонтальный, то возвращаемое значение соответствует значению K.
int ntape(int line,int dir)
{
	int m;
	switch (dir)
	{
	case 0:
		m=line;
		break;
	case 1:
		m=line-1;
		break;
	case -1:
		m=line;
	}
return m;
}

//======================================================================
// powers of j(j+1) and k to be root squared in the four tapes of diagram;
// values of <k|Jz|k> for five lines.
//----------------------------------------------------------------------
void Fill_tapes(Tape& pjk_tape, Tape& pJz)
{
	j_k_powers J,K2,K1,C;
	Term Tm;
	Tape A;

// j(j+1)
	J.jp=1;
	J.kp=0;
	J.val=1;
// K^2
	K2.jp=0;
	K2.kp=2;
	K2.val=-1;
// K^1
	K1.jp=0;
	K1.kp=1;
	K1.val=-1;
// Const
	C.jp=0;
	C.kp=0;
	C.val=0;

// Template for one term to be root squared (for one up elementary path)
	Tm.push_back(J);
	Tm.push_back(K2);
	Tm.push_back(K1);
	Tm.push_back(C);

// <k-1||k>
	A.push_back(Tm);
	A[0][2].val=1;
	A[0][3].val=0;
// <k||k+1>
	A.push_back(Tm);
// <k+1||k+2>
	A.push_back(Tm);
	A[2][2].val=-3;
	A[2][3].val=-2;
// <k+2||k+3>
	A.push_back(Tm);
	A[3][2].val=-5;
	A[3][3].val=-6;
// <k+2||k+3>
	A.push_back(Tm);
	A[4][2].val=-7;
	A[4][3].val=-12;

	pjk_tape=A;

// K^1
	K1.jp=0;
	K1.kp=1;
	K1.val=1;
	Tm.clear();

// Template for <k|Jz|k>)
	Tm.push_back(K1);
	Tm.push_back(C);
	A.clear();
	for(int k=0; k<5; k++)
	{
		A.push_back(Tm);
		A[k][1].val=k-2;
	}

	pJz=A;

}

//======================================================================
//----------------------------------------------------------------------
void print_Tape(Tape pJz, char* s)
{
	ff << s << endl;
	for(int k=0; k<pJz.size(); k++)
	{
		ff << "Term " << k << endl;
		for(int i=0; i<pJz[k].size(); i++)
		{
			ff.width(4);
			ff << pJz[k][i].jp;
			ff.width(4);
			ff << pJz[k][i].kp;
			ff.width(10);
			ff << pJz[k][i].val << endl;
		}
	}
}

//======================================================================
//----------------------------------------------------------------------
void print_Term(Term T, char* s)
{
	ff << s << endl;
	for(int i=0; i<T.size(); i++)
	{
		ff.width(4);
		ff << T[i].jp;
		ff.width(4);
		ff << T[i].kp;
		ff.width(10);
		ff << T[i].val << endl;
	}
}
//======================================================================
//----------------------------------------------------------------------
Term Terms_multiplication(Term A,Term B)
{
	Term C;
	j_k_powers w;

	C.clear();
	for(int k=0; k<A.size(); k++)
		for(int i=0; i<B.size(); i++)
		{
			w.jp=A[k].jp+B[i].jp;
			w.kp=A[k].kp+B[i].kp;
			w.val=A[k].val*B[i].val;
			C.push_back(w);
		}
	return C;
}

//======================================================================
// Summation of the two terms. Each power of j and k in the term 
// must be unique.
//----------------------------------------------------------------------
Term Terms_summation(Term A,Term B)
{
	Term C;
	j_k_powers w;

//print_Term(A,"A");
//print_Term(B,"B");

	C.clear();
// Common terms from A and B.
	for(int k=0; k<A.size(); k++)
		for(int i=0; i<B.size(); i++)
		{
			if((A[k].jp==B[i].jp) && (A[k].kp==B[i].kp))
			{
				w=A[k];
				w.val+=B[i].val;
				C.push_back(w);
			}
		}
// Unique terms from A
	for(int i=0; i<A.size(); i++)
	{
		int y=0;
		for(int k=0; k<B.size(); k++)
			if((B[k].jp==A[i].jp) && (B[k].kp==A[i].kp))
				y=1;
		if(y==0)
			C.push_back(A[i]);
	}
// Unique terms from B
	for(i=0; i<B.size(); i++)
	{
		int y=0;
		for(int k=0; k<A.size(); k++)
			if((A[k].jp==B[i].jp) && (A[k].kp==B[i].kp))
				y=1;
		if(y==0)
			C.push_back(B[i]);
	}
	sort(C.begin(),C.end(),compj_k_powers);
	EqualAdds(C);
	return C;
}

//======================================================================
// Логическая функция для сортировки вектора Term.
//----------------------------------------------------------------------
bool compj_k_powers (j_k_powers A, j_k_powers B)
{
	return (wmn*A.jp+A.kp > wmn*B.jp+B.kp);
}

//======================================================================
// Сбор подобных членов.
//----------------------------------------------------------------------
void EqualAdds(Term& A)
{
	Term B=A,C;
	j_k_powers w;

	sort(B.begin(),B.end(),compj_k_powers);
//print_Term(B,"Sorted path");
	C.clear();
	int i=0, k=0;
	while(i<B.size())
	{
		w=B[i];
//ff << i << endl;
		C.push_back(w);
		i++;
		while((i<B.size()) && (w==B[i]))
		{
			C[k].val+=B[i].val;
//ff << i << " " << k << endl;
			i++;
		}
		k++;
	}
//	print_Term(C,"Ultimate one path");
// Deleting terms with zero coefficients/
	k=0;
	B.clear();
	for(i=0; i<C.size(); i++)
		if(C[i].val != 0)
			B.push_back(C[i]);
	A=B;
}

bool j_k_powers::operator==(const j_k_powers& A) const
{
	return(this->jp == A.jp && this->kp == A.kp);
}

//**************************************
//     
// Name: binomial coefficiant recursion
// Description:calculate binomial coeffi
//     ciants recursively
// By: Anthony Bernardino
//
// Inputs:into 2 integers, 1st integer m
//     ust be larger than the 2nd integer.
//
// Returns:returns binomial coefficiant
//
// Assumes:program is recursive, so be c
//     areful.
//
//This code is copyrighted and has// limited warranties.Please see http://
//     www.Planet-Source-Code.com/vb/scripts/Sh
//     owCode.asp?txtCodeId=6237&lngWId=3//for details.
//     

int binom (int n, int i)
{
    int n1;
    int n2;
    if (( i ==0 ) || (n== i))
    {
        return 1;
    }
    else
    {
        n1 = binom ( n - 1, i);
        n2 = binom ( n - 1, i -1);
        return n1 + n2;
    }
}

char fs(int i)
{
	switch(i)
	{
	case(0):
	return 'x';
	break;

	case(1):
	return 'y';
	break;

	case(2):
	return 'z';
	break;

	}
}

Jway.cpp Implementation of the Jway class.

// Jway.cpp: implementation of the Jway class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "Jway.h"
#include <iostream>

using namespace std;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Jway::Jway()
{
}

Jway::~Jway()
{
}

void Jway::Init(rlist r)
{
	switch (r)
	{
	case x:
		nways=2;
		way.push_back(-1); //  <k||k-1>
		way.push_back(1);  //  <k||k+1>
		break;

	case y:
		nways=2;
		way.push_back(-1); //  <k||k-1>
		way.push_back(1);  //  <k||k+1>
		break;

	case z:
		nways=1;
		way.push_back(0);  //  <k||k>
		break;

	default:
		cout << "The argument in Jway::Init is false" << endl;
		exit(1);
	}
}


2. Вывод аналитических выражений, связывающих компоненты обобщенного оператора углового момента (неприводимый тензорный оператор) и декартовые компоненты этого же оператора

Эффективные вращательные гамильтонианы многоатомных молекул обычно выражаются через степенные ряды декартовых компонент оператора углового момента. К примеру, следующим образом (формула Уотсона):

(1)

Здесь h - параметры, J - оператор углового момента.

Если молекула обладает какой-либо симметрией, то часть членов в сумме (1) пропадает. Это связано с тем, что H_eff должен оставаться неизменным при всех операциях симметрии. Операциями симметрии являются повороты и отражения системы координат, которые совмещают эквивалентные атомы (то есть атомы одного сорта). Выявить, какие члены остаются, а какие пропадают в гамильтониане, записанном в виде (1), довольно сложно, так как декартовые компоненты при поворотах выражаются друг через друга и к тому же они не коммутируют. Сделать это намного проще, если записывать гамильтониан через обобщенные операторы углового момента (ООМ), которые преобразуются по неприводимым представлениям группы вращений O^+(3). Такие операторы еще называются неприводимыми тензорными операторами.

Определение обобщенного оператора углового момента.

Неприводимый тензорный оператор (НТО) характеризуется рангом L, который может быть любым целым положительным числом начиная с 0. НТО ранга L имеет 2L+1 компоненту. Каждая компонента имеет свой номер M, который может принимать любые целые значения в диапазоне от -L до L. Обозначение любой компоненты выглядит следующим образом: T_M^L. Основное соотношение, которому должны удовлетворять компоненты НТО, следующее:

Здесь J_mu циклические компоненты оператора углового момента, ccg коэффициенты Клебша - Гордона.

НТО могут быть любой природы. Если они сконструированы из компонент обычного углового момента, то тогда они называются обобщенными операторами углового момента (ООМ).

Эффективный вращательный гамильтониан, выраженный через ООМ, будет выглядеть следующим образом:

(2)

Гамильтониан в такой форме полностью эквивалентен (1), но значительно лучше приспособлен для учета симметрии молекул. В практических задачах необходимо уметь преобразовывать гамильтониан из формы (1) в форму (2) и наоборот. Для этого нужны аналитические выражения, связывающие ООМ с декартовыми компонентами углового момента.

Целью представленной ниже программы является вывод таких выражений.

Исходным выражением для решения поставленной задачи является рекуррентное соотношение между тремя компонентами ООМ с одинаковым номером и разными рангами:

(3)

Алгоритм построения выражения для произвольного T_M^L через операторы J^2, J_x, J_y, J_z состоит в последовательном применении формулы (3) с возрастающими значениями L. При этом базовые операторы T_M^L с минимально возможным L, когда L=M, определяются формулой

Программа, используя рекурентное соотношение (3), ищет T_M^L в наиболее общем виде, который имеет вид

(4)

Здесь - целые числа.

При своей работе при перемножении двух операторных многочленов программа собирает подобные члены, ищет общие множители в числителе и знаменателе и сокращает их, извлекает корни из множителей, которые представляют собой квадраты целых чисел и т.д.

Программа состоит из двух файлов, заголовочного types.h и исполняемого main.cpp.

types.h Определяет составные типы данных, используемые в данной задаче, и объявляет прототипы функций.

#if !defined(types_h)
#define types_h

#include <vector>
#include <list>
#include <string>

using namespace std;

//typedef struct
class j_k_powers
{
public:
	double val; // factor;
	int jp;  // power of J^2;
	int kp;  // power of J_z;
    bool operator==(const j_k_powers& ) const;
};

typedef vector<j_k_powers> Term; /* regular expression for Tlm */

class Tlm {
public:
	Term clm;
	list<int> a; /* Простые множители числителя под корнем */
	list<int> b; /* Простые множители знаменателя под корнем */
	/* list используется, чтобы было легче упрощать a/b */
};

typedef list<int> lstarr;


void PrimeFactors(int ,vector<int>& );
Term Terms_multiplication(Term ,Term );
void EqualAdds(Term& );
bool compj_k_powers (j_k_powers , j_k_powers );
void division(list<int>& , list<int>& );
bool comp_int (int , int );
list<int>  lists_devide(list<int>& , list<int>& );
list<int> lists_mult(list<int> , vector<int> );
void print_vect(const string,list<int> );
void clear(Tlm& );
string IntToString (int );
int SqrtInt(list<int> );
int fullnumber(list<int>);

void list_erase(void);

#endif //types_h

main.cpp Основная и вспомогательные фйнкции.

// main.cpp : Defines the entry point for the console application.
// Program to evaluate expressions for irreducible tensor operators.
// Next version after tlm.cpp. Using rearanged formulas.

#include "types.h"
#include <fstream>
#include <algorithm>
#include <string>
#include <sstream>
#define wmn 1000  // to sort vector

using namespace std;
ofstream ff;

int main(int argc, char* argv[])
{
	ff.open("out.txt",ios::out);
	int M = 4;  /* Значение M цепочки TLM */
	int QuantityOfTlm = 12;

	ff << "Рассчитываются TLM для M=" << M << endl;

	vector<Tlm> TlmArray; /* Массив для всех найденных Tlm. Первую позицию занимает TMM */
	TlmArray.reserve(QuantityOfTlm);

	string s1 = IntToString(M);
	string s2;

	int Lmax = M + QuantityOfTlm;

// Рабочие переменные
	Tlm tmw;
	j_k_powers w;
	vector<int> f; /* prime factors */

	f.resize(40);

// == T^(M)_(M) =========================
// T^(M)_(M) has just one term. T^(M)_(M) = 1*(T^(M)_(M)).
	ff << endl << "=========== " << "T^"+s1+"_"+s1 << " ===========" << endl;
	clear(tmw);
	w.val = 1;
	w.jp = 0;
	w.kp = 0;
	tmw.clm.push_back(w);
	tmw.a.push_back(1);
	tmw.b.push_back(1);

	TlmArray.push_back(tmw);

	ff << "val " << "jp " << "kp" << endl;
	ff << w.val << " " << w.jp << " " << w.kp << endl;

// --------------------------------------

// == T^(M+1)_(M) =======================
	s2 = IntToString(M+1);
	ff << endl << "=========== " << "T^"+s2+"_"+s1 << " ===========" << endl;
	clear(tmw);

	PrimeFactors(2*M+1,f); // Раскладывает на простые множители.
// Вставляем обязательно 1, для того, чтобы при сокращении множителей
// не получить нулевого результата.
// Числитель
	tmw.a.push_back(1);
	for(int i=0; i < f.size(); i++)
	{
		tmw.a.push_back(f[i]);
	}
// Знаменатель
	tmw.b.push_back(1);
	tmw.b.push_back(2);
	tmw.b.push_back(2);

	ff << "val " << "jp " << "kp" << endl;

	w.val = M;
	w.jp = 0;
	w.kp = 0;
	if(M != 0){
		tmw.clm.push_back(w);
		ff << w.val << " " << w.jp << " " << w.kp << endl;
	}

	w.val = 2;
	w.jp = 0;
	w.kp = 1;
	tmw.clm.push_back(w);
	ff << w.val << " " << w.jp << " " << w.kp << endl;

	TlmArray.push_back(tmw);
// ------------------------------------------

// == all Tlms с использованием рекурентной формулы ==
// Соответсвия: i=0 -> T^(M)_(M), i=1 -> T^(M+1)_(M), i=2 -> T^(M+2)_(M),....
	Tlm tmw1, tmw2; /* две составляющие нового Tlm */
	Term r, g1, g2;
	for(i = 2; i < QuantityOfTlm; i++)
	{
		clear(tmw1);
		clear(tmw2);
		r.clear();

		int L=M+i;  /* очередное L, для которого строим выражение */
		s2 = IntToString(L);
		ff << endl << "=========== " << "T^"+s2+"_"+s1 << " ===========" << endl;
		int s = i-1; /* номер последнего вычисленного Tlm */

// Включаем в выражение член с T^{L-1}_{M} -------------------

		w.val = M;
		w.jp = 0;
		w.kp = 0;
		r.push_back(w);
		w.val = 2;
		w.jp = 0;
		w.kp = 1;
		r.push_back(w);
/*
* Умножение P_{L-1}*(2J_z+M) и приведение подобных.
*/
/* g1 содержит полином от степеней J^2 и J_z в первой части формулы,
то есть той, которая содержит T^{L-1}_{M}  */
		g1 = Terms_multiplication(TlmArray[s].clm,r);
		EqualAdds(g1);

/* tmw1.a, tmw1.b - числитель и знаменатель подкоренного выражения в первой части формулы,
то есть той, которая содержит T^{L-1}_{M}  */
// Знаменатель добавляемого подкоренного выражения
		PrimeFactors((L+M)*(L-M)*16,f);
		tmw1.b = lists_mult(TlmArray[s].b, f);
// Числитель добавляемого подкоренного выражения
		PrimeFactors((2*L-1)*(2*L-1)*4,f);
		tmw1.a = lists_mult(TlmArray[s].a, f);
		int j1;

// Член в выражении с T^{L-2}_{M} -------------------

		r.clear();
		w.val = L*(L-2);
		w.jp = 0;
		w.kp = 0;
		r.push_back(w);
		w.val = -4;
		w.jp = 1;
		w.kp = 0;
		r.push_back(w);
/* g2 содержит полином от степеней J^2 и J_z во второй части формулы,
то есть той, которая содержит T^{L-2}_{M}  */
		g2 = Terms_multiplication(TlmArray[s-1].clm,r);
		EqualAdds(g2);

/* tmw2.a, tmw2.b - числитель и знаменатель подкоренного выражения во второй части формулы,
то есть той, которая содержит T^{L-2}_{M}  */
// Знаменатель добавляемого подкоренного выражения
		PrimeFactors((L+M)*(L-M)*16,f);
		tmw2.b = lists_mult(TlmArray[s-1].b, f);
		PrimeFactors((L-1+M)*(L-1-M),f);
		tmw2.a = lists_mult(TlmArray[s-1].a, f);
		
	list<int> bl = lists_devide(tmw1.b, tmw2.b);
/* bl - общий множитель подкоренного выражения.
В tmw1.b, tmw2.b остаются индивидуальные множители */

/* Анализируем и приводим к общему знаменателю, используя bl, tmw1.b и tmw2.b*/
		for(list<int>::iterator k = tmw1.b.begin(); k != tmw1.b.end(); k++)
			bl.push_back(*k);
		for(k = tmw2.b.begin(); k != tmw2.b.end(); k++)
			if(*k != 1)  
				bl.push_back(*k);
		bl.sort();
/* bl содержит общий знаменатель */

/* Числитель подкоренного выражения при T^{L-1}_{M}*/
		for(k = tmw2.b.begin(); k != tmw2.b.end(); k++)
			if(*k != 1) 
				tmw1.a.push_back(*k);
		tmw1.a.sort();
/* Числитель подкоренного выражения при T^{L-2}_{M}*/
		for(k = tmw1.b.begin(); k != tmw1.b.end(); k++)
			if(*k != 1) 
				tmw2.a.push_back(*k);
		tmw2.a.sort();

		list<int> a1=tmw1.a, a2=tmw2.a;
// bc - общий множитель числителей под корнем, который нужно вынести за скобки
// и попытаться сократить.
		list<int> bc = lists_devide(a1, a2);
		int sqa1 = SqrtInt(a1);
		int sqa2 = SqrtInt(a2);
		if(sqa1*sqa1 != fullnumber(a1))
		{
			ff << "Какие-то проблемы. 1" << endl;
			exit(1);
		}
		if(sqa2*sqa2 != fullnumber(a2))
		{
			ff << "Какие-то проблемы. 2" << endl;
			exit(2);
		}
//Корректировка левого полинома
		for(int j=0; j < g1.size(); j++)
			g1[j].val*=sqa1;

//Корректировка правого полинома
//		ff << "g2s " << g2.size() << endl;
		for(j=0; j < g2.size(); j++)
			g2[j].val*=sqa2;

// Сложение левого и правого полиномов с результатом в g1.
		for(j=0; j < g2.size(); j++)
			g1.push_back(g2[j]);

// Сбор подобных членов в полиноме.
		EqualAdds(g1);

// Определение общего множителя в полиноме
		vector<lstarr> mutual;
		mutual.clear();
		list<int> wlist;
		for(j=0; j < g1.size(); j++)
		{
			f.clear();
			wlist.clear();
			PrimeFactors(g1[j].val,f);
			for(vector<int>::iterator k = f.begin(); k != f.end(); k++)
				wlist.push_back(*k);
			mutual.push_back(wlist);

		}

		wlist=mutual[0];
		for(j=1; j < mutual.size(); j++)
			wlist=lists_devide(wlist, mutual[j]);
		
		int icf = fullnumber(wlist);
		tmw.clm = g1;
		ff << "Результирующий полином" << endl;
		for(j=0; j < tmw.clm.size(); j++)
		{
			tmw.clm[j].val/=icf;
			ff << tmw.clm[j].val << " " << tmw.clm[j].jp << " " << tmw.clm[j].kp << endl;
		}

		tmw.a = bc;
		tmw.b = bl;
// Переносим общий множитель полинома в числитель множителя всего выражения под корень
		for(list<int>::iterator j2=wlist.begin(); j2 != wlist.end(); j2++)
		{
			if(*j2 != 1)
			{
				tmw.a.push_back(*j2);
				tmw.a.push_back(*j2);
			}
		}
		lists_devide(tmw.a, tmw.b);
	print_vect("Числитель подкоренного выражение",tmw.a);
	print_vect("Знаменатель подкоренного выражение",tmw.b);
		TlmArray.push_back(tmw);
	}

	return 0;
}

//======================================================================
/* ===== Разложение числа на простые множители ===== */
//----------------------------------------------------------------------
void PrimeFactors(int s,vector<int>& f)
{
	f.clear();
	int t, i=2;
// Число, раскладываемое на простые множители должно быть положительным.
	t = abs(s);
	while(i<=t)
	{
		if(t%i==0)
		{
			f.push_back(i);
			t=t/i;
		}
		else
			i=i+1;		
	}
}

//======================================================================
// Умножение двух полиномов от степеней J^2 и J_z 
//----------------------------------------------------------------------
Term Terms_multiplication(Term A,Term B)
{
	Term C;
	j_k_powers w;

	C.clear();
	for(int k=0; k<A.size(); k++)
		for(int i=0; i<B.size(); i++)
		{
			w.jp=A[k].jp+B[i].jp;
			w.kp=A[k].kp+B[i].kp;
			w.val=A[k].val*B[i].val;
			C.push_back(w);
		}
	return C;
}

//======================================================================
// Сбор подобных членов
//----------------------------------------------------------------------
void EqualAdds(Term& A)
{
	Term B=A,C;
	j_k_powers w;

	std::sort(B.begin(),B.end(),compj_k_powers);
//print_Term(B,"Sorted path");
	C.clear();
	int i=0, k=0;
	while(i<B.size())
	{
		w=B[i];
//ff << i << endl;
		C.push_back(w);
		i++;
		while((i<B.size()) && (w==B[i]))
		{
			C[k].val+=B[i].val;
//ff << i << " " << k << endl;
			i++;
		}
		k++;
	}
//	print_Term(C,"Ultimate one path");
// Deleting terms with zero coefficients/
	k=0;
	B.clear();
	for(i=0; i<C.size(); i++)
		if(C[i].val != 0)
			B.push_back(C[i]);
	A=B;
}

//======================================================================
// Логическая функция сравнения для сортировки вектора Term.
//----------------------------------------------------------------------
bool compj_k_powers (const j_k_powers A, const j_k_powers B)
{
	int a = wmn*A.jp+A.kp;
	int b = wmn*B.jp+B.kp;
	return (a > b);
//	return (a >= b);
//	return (wmn*A.jp+A.kp >= wmn*B.jp+B.kp);
}

//======================================================================
// Оператор равенства степеней.
//----------------------------------------------------------------------
bool j_k_powers::operator==(const j_k_powers& A) const
{
	return(this->jp == A.jp && this->kp == A.kp);
}

//======================================================================
// Деление двух целых чисел представленных в виде произведения простых чисел.
//----------------------------------------------------------------------
void division(std::list<int>& a, std::list<int>& b)
{
//	sort(a.begin(),a.end(),comp_int);
//	sort(b.begin(),b.end(),comp_int);
	a.sort();
	b.sort();
}

//======================================================================
// Логическая функция сравнения для сортировки вектора std::vector<int>.
//----------------------------------------------------------------------
bool comp_int (int A, int B)
{
	return (A >= B);
}

//======================================================================
// Исключение общих простых множителей из двух списков простых множителей.
// a - числитель.
// b - знаменатель.
// Общие множители - возвращаемое значение.
//------------------------------------------------------------------
list<int> lists_devide(list<int>& a, list<int>& b)
{
    list<int> o;

	for( list<int>::iterator k = b.begin(); k != b.end();)
	{
		int c = 0;
	    for( list<int>::iterator i = a.begin(); i != a.end();)
		{
			if(*k == *i && *k != 1)
			{
				i=a.erase(i);
				c = 1;
				break;
			}
			else
				i++;
		}
		if(c == 1)
		{
			o.push_back(*k);
			k=b.erase(k);
/*
			ff<< "a" << endl;
			for( i = a.begin(); i != a.end();i++)
				ff<<*i<<" ";
			ff<<endl;
			ff<< "b" << endl;
			for( i = b.begin(); i != b.end();i++)
				ff<<*i<<" ";
			ff<<endl;
*/
		}
		else
			k++;
	}
	return o;
}

//======================================================================
// Умножение (объединение) двух строк простых множителей.
//------------------------------------------------------------------
list<int> lists_mult(list<int> a, vector<int> b)
{
	list<int> c;
	c.clear();
	c = a;
	for(vector<int>::iterator i = b.begin(); i != b.end(); i++)
		c.push_back(*i);
	c.sort();
	return c;
}

//======================================================================
// Перемножение всех членов list<int>.
//------------------------------------------------------------------
int fullnumber(list<int> a)
{
	int m=1;
    for( list<int>::iterator k = a.begin(); k != a.end();k++)
		m =m* (*k);
	return m;
}

/* ======================================================================
* Извлечение корня из целого числа представленного в виде произведения
* простых множителей и являющегося квадратом целого числа.
* В функции выполняется, что число действительно есть квадрат целого числаю
* Проверка осуществляется путем определения количества простых множителей.
------------------------------------------------------------------ */
int SqrtInt(list<int> f)
{
	list<int> fin = f;
	fin.sort();
	int c=0, o=1;
	for(list<int>::iterator i=fin.begin(); i != fin.end();)
	{
		int m;
		int k;  // количество одинаковых
		m=*i; k=0;
		while(m == *i && i != fin.end())
		{
			k++; i++;
		}
		if(k%2 != 0 && m != 1)
		{
			ff << "Корень извлечь нельзя" << endl;
			return -1;
		}
		k=k/2;
		for(int l=0; l<k; l++)
			o*=m;
	}
	return o;
}


//======================================================================
// Очистка полей Tlm.
//------------------------------------------------------------------
void clear(Tlm& tmw)
{
	tmw.a.clear();
	tmw.b.clear();
	tmw.clm.clear();
}


//======================================================================
// Преобразование int to string.
//------------------------------------------------------------------
string IntToString (int a)
{
    string str;
    ostringstream temp;
    temp<<a;
    return temp.str();
}
//======================================================================
// Вспомогательный код для отработки разных алгоритмов.
//------------------------------------------------------------------

void list_erase(void)
{
	list<int> x;
    for ( int i =0;i<10; i++)
        x.push_back(i);

    for( list<int>::iterator k = x.begin(); k != x.end();k++)
        ff<<*k<<" ";
    ff<<endl;

	for( k = x.begin(); k != x.end();)
	  if( !((*k)%2) )
		k=x.erase(k);
	  else
		++k;

    for( k = x.begin(); k != x.end();k++)
        ff<<*k<<" ";
    ff<<endl;
}

//======================================================================
// Вывод в поток ff компонент list<int> f.
// h = 1,2 - два разных формата вывода.
//------------------------------------------------------------------
void print_vect(const string inf,list<int> f)
{
	int h=1;
//	ff << inf << endl;
	ff << inf+" = " << fullnumber(f) << "  ->  ";
	if(h==1)
	{
		for(list<int>::iterator i=f.begin(); i != f.end(); i++)
		{
			ff << *i << " ";
		}
		ff << endl;
	}
	h=2;
	if(h==2)
	{
// Извлечение корня, а затем вывод в поток	
		list<int> o, fin = f;
		fin.sort();
		int c=0;
		for(list<int>::iterator i=fin.begin(); i != fin.end();)
		{
			int m=*i;
			int k=0;  // количество одинаковых
			while(m == *i && i != fin.end())
			{
				k++; i++;
			}
			ff << m << ", ко-во =" << k << endl;

		}
	}
}

Фрагмент распечатки результатов работы программы.

=========== T^5_3 ===========

Результирующий полином

val jp kp
-1 1 0
9 0 2
27 0 1
24 0 0

Числитель подкоренного выражения = 7 -> 7

7, к-во =1

Знаменатель подкоренного выражения = 16 -> 1 2 2 2 2

1, к-во =1

2, к-во =4

Эти же результаты, записанные в виде формулы, имеют вид: